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ABSTRACT 

A  summary  is  given  of  significant  results  obtained  during  the  research 
effort  corresponding  to  the  period  1  September  1968  -  30  September  1969. 
Particular  emphasis  was  given  to  the  construction  of  a  phase-path  doppler 
sounder  array.  Signals  detected  by  the  microbarograph  array  on  the  occa¬ 
sion  of  the  Chinese  nuclear  test  of  29  September  1969  are  compared  with 
signals  from  previous  tests.  There  appears  to  be  a  seasonal  effect  on  the 
propagation  and/or  bandwidth  of  these  signals.  Preliminary  work  on  ionospheric 
motion  background  using  phase-path  sounder  data  indicates  that  the  background 
activity  can  be  broadly  separated  into  two  period  ranges:  periods  longer  and 
shorter  than  5  minutes.  The  shorter  periods  appear  to  have  a  hydromagnetic 
origin  while  the  longer  periods  are  probably  due  to  internal  gravity  wave 
activity.  A  study  of  the  spatial  coherency  of  the  signal  and  noise  field  for 
the  array  of  microbarographs  indicated  that  acoustic  signals  have  coherencies 
well  above  the  background  noise  at  separations  greater  than  3  or  4  km.  The 
signal  coherency  remains  high  up  to  separations  of  50  km  and  transverse 
signal  coherencies  are  larger  than  parallel  signal  coherencies  as  would  be 
expected  from  long-crested  waves.  A  theoretical  investigation  of  the  inter¬ 
action  between  instability  waves  and  internal  gravity  waves  in  the  atmosphere 
showed  that  internal  gravity  waves  are  trapped  within  a  layer  where  the  shear 
flow  velocity  is  greater  than  the  speed  of  sound  and  that  the  phase  velocity 
of  the  internal  gravity  waves  is  approximately  equal  to  the  maximum  velocity 
of  the  shear  flow.  As  part  of  the  theoretical  aspect  of  signal  processing 
two  array  processing  techniques  have  been  developed. 
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1.  INTRODUCTION 

In  the  course  of  this  year’s  research  effort,  the  authors  have  had 
occasion  to  make  a  critical  review;  of  the  course  of  the  project,  its  history, 
previous  results,  and  what  they  felt  to  be  its  future  direction.  This  apprai¬ 
sal  has  led  them  in  some  cases  to  redirect  or  shift  the  orientation  and 
emphasis  of  the  research  in  the  light  of  previous  and  ne f  results. 

i 

It  is  pertinent  to  review  the  prior  history  of  the  infrasonics 
program  since  most  of  its  i  litial  principals  have  departed  the  scene.  The 
project  started  in  late  1965  at  the  Hudson  Laboratories  of  Columbia  University, 
being  funded  by  ARPA  and  administered  by  the  Office  of  Naval  Research.  It  is 
our  understanding  that  the  underlying  purpose  of  the  project  was  to  look  for 
the  existence  of  long  period  (greater  than  10  minutes)  gravity  waves  in  the 
atmosphere  and  to  study  their  possible  usefulness  in  the  detection  and  diagnosis 
of  atmospheric  nuclear  detonations.  The  initial  phases  were  devoted  to  the 
design  and  construction  of  suitable  sensors  which  were  installed  in  remote 
locations  about  the  middle  of  1966.  A  paper- tape  digital  data  acquistion 
system  was  installed  initially  and  computer  programs  for  data  processing  were 
written.  The  Chinese  atmospheric  test  of  June  1967  gave  some  encouragement  to 
the  validity  of  the  concept  by  a  definitive  recording.  The  data  acquisition 
system  was  later  improved  so  that  the  output  could  be  written  on  magnetic  tape 
in  a  form  compatible  with  the  subsequent  data  processing.  Prior  to  the  French 
series  of  atmospheric  tests  in  the  summer  of  1968,  the  microbarograph  array  had 
grown  to  some  12  instruments  and  had  an  aperture  of  about  250  kilometers. 

Later  a  single  phase-path  doppler  ionosonde  and  a  three-component  magnetometer 
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array  were  installed.  The  magnitude  of  the  arrays  brought  with  it  concomitant 
financial  outlays  in  maintenance,  personnel,  equipment  and  telephone  rental. 

Also  at  this  time  the  project  used  an  IBM-1800  computer,  having  had  an  1130 
computer  with  limited  capacity. 

In  April  1968  Columbia  University  announced  the  impending  closure  of 
Hudson  Laboratories.  At  that  time  arrangement  were  made  to  transfer  the  pro¬ 
ject  to  Teledyne  Isotopes  as  of  September  1968.  After  a  period  of  shutdown 
during  the  transfer,  extended  because  of  the  prolonged  AT&T  strike  of  1968, 
the  system  again  became  operational  at  its  new  location  in  Westwood .  An  IBM 
360/30  modified  by  software  to  operate  in  the  360/44  mode  for  scientific  work 
was  installed  at  Westwood  to  replace  the  IBM-1800  computer.  Both  the  hardware 
and  software  have  proved  entirely  satisfactory  for  all  requirements  of  the 
project,  and  no  modification  in  either  is  being  sought.  The  Hudson  Laboratories 
computer  programs  were  translated  to  the  IBM-360  language  and  some  new  programs 
have  been  added  to  the  project’s  computer  library. 

During  the  period  from  mid-1966  to  the  present,  theoretical  studies 
and  data  analyses  have  been  carried  on  concurrently  with  the  data  gathering. 

The  results  of  these  analyses  and  supporting  data  have  beta  presented  in  various 
technical  and  quarterly  progress  reports  and  may  be  summarized  as  follows. 

Definite  observations  of  long  period  gravity  waves  have  been  made  in 
conjunction  with  atmospheric  nuclear  detonations. 

The  pressure  signals  themselves  are  rather  weak.  The  signal-to-noise 
ratio  at  any  particular  station  is  generally  less  than  one.  It  follows,  there¬ 
fore,  that  the  signals  can  only  be  extracted  from  the  noise  by  judicious  filtering. 
The  data  from  any  station  where  the  noise  power  is  considerably  in  excess  of  the 
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mean  noise  power,  must  be  eliminated.  The  data  must  also  be  filtered  with  a 
fairly  narrow  band-pass  filter.  Further,  the  array  must  be  beam  steered,  i.e. 
by  velocity  and  directional  filtering,  and  a  knowledge  of  Tq  is  essential  to 
identification . 

Theoretical  studies  have  shown  that  a  large  part  of  the  energy  of 
the  long  period  gravity  waves,  generated  by  nuclear  explosions,  will  be  trapped 
near  the  "top”  of  the  atmosphere.  Particle  motion  should  therefore  be  large 
at  ionospheric  heights  and  the  corresponding  pressure  fluctuations  should  be 
small  at  ground  level.  Such  behavior  appears  to  be  verified  by  the  observa¬ 
tion  of  a  signal  using  a  phase-path  doppler  sounder  which  we  believe  to  be  a 
gravity  wave  generated  by  the  French  nuclear  detonation  of  24  August  1968. 

The  signal- to-noise  ratio  for  the  observed  event- related  signal  is  about  2. 
Since  only  one  ionosonde  was  in  operation,  no  information  could  be  obtained 
concerning  the  phase  velocity  and  direction  of  travel  of  this  wave. 

On  the  basis  of  such  results  we  were  forced  to  conclude  that  micro¬ 
barographs  do  not  appear  to  be  the  best  sensors  with  which  to  detect  long 
period  gravity  waves  for  the  program's  expressed  purpose  of  detection  and 
diagnosis.  Furthermore,  a  comparison  of  the  theoretical  response  of  the 
Teledyne-Isotopes  array  with  that  of  a  number  of  other  arrays  such  as  10, 

15,  and  21-element  isometric  arrays  and  the  LAMA  array,  showed  that  not  much 
improvement  could  be  expected  by  adding  new  instruments  or  changing  the  geom¬ 
etry  of  the  present  array.  It  would  appear  rather  that  the  phase-path  doppler 
ionosonde  is  perhaps  a  more  effective  sensor.  On  this  premise  we  considered 
concentrating  further  effort  on  the  phase-path  sounder. 
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The  phase-path  doppler  sounder  array  is  only  in  its  initial  experi¬ 
mental  phase,  and  much  remains  to  be  explored  in  relation  to  optimum  transmission 
frequencies,  noise  and  signal  levels,  and  coherence  over  the  array.  Placing 
greater  emphasis  on  the  doppler  array  has  introduced  a  new  class  of  ionospheric 
phenomena  not  relevant  to  ground  level  pressure  measurements.  These  include, 
for  instance,  noise  associated  with  solar  flares,  geomagnetic  storms,  hydro- 
magnetic  modulation  of  the  ionospheric  plasma,  interaction  between  neutral 
and  ionized  components  of  the  upper  atmosphere. 

On  the  occasion  of  the  Chinese  nuclear  test  of  29  September  1969, 
signals  coming  from  the  direction  of  the  test  site  were  recorded  in  the  micro¬ 
barograph  array.  The  analyses  of  these  recordings  were  carried  out  in  a 
similar  way  as  the  signals  from  previous  tests  were  analyzed  at  Hudson  Labora¬ 
tories,  so  a  direct  comparison  between  previous  results  and  results  from  this 
test  can  be  made. 

As  part  of  the  phase-path  variation  studies,  preliminary  work  has 
been  done  on  ionospheric  background  phenomena.  This  work  has  been  directed 
mainly  to  the  determination  of  the  type  of  driving  mechanism  characterizing 
various  period  ranges. 

Experimental  results  at  Hudson  Laboratories  indicated  that  internal 
gravity  waves  are  associated  with  the  jet  stream.  It  has  been  suggested  that 
these  internal  waves  are  generated  by  a  non-linear  interaction  of  instability 
waves.  It  was  therefore  desirable  to  carry  out  a  theoretical  study  of  the 
interaction  of  instability  waves  and  internal  gravity  waves.  A  model  which 
neglects  compressibility  was  used  to  find  solutions  for  this  problem. 
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Finally  the  important  task  concerned  with  the  recognition  of  signals 
from  the  noise  has  been  pursued  by  Dr.  M.  Hinich  of  the  Carnegie  Mellon  Insti¬ 
tute  and  Dr.  E.  Posmentier  of  we w  York  University,  resulting  in  two  array  pro¬ 
cessing  techniques. 
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2.  EXPERIMENTAL  ATMf  °HERIC  PHYSICS 


2.1  Teledyne  Isotopes  Array 
2.1.1  Array  Geography 

Data  for  atmospheric-wave  propagation  studies  is  collected  by  an 
array  of  sensors,  consisting  of  phase-path  doppler  sounders,  microbarographs 
and  magnetometers,  located  at  the  various  stations  shovm  in  Figure  2.1  .  The 
geographic  coordinates  of  these  stations  are  listed  in  Table  2.1. 

Phase-path  doppler  sounder  array 

A  continuous  wave  4-element  doppler  sounder  array  has  been  in 
operation  during  the  last  months  of  the  contract  year.  The  transmitter  is 
situated  at  Sterling  Forest.  It  should  be  noted  that  this  is  the  only  location 
for  which  it  has  been  possible  thus  far  to  obtain  a  transmission  license  from 
the  F.C.C.  The  array  is  licensed  to  operate  at  frequencies  of  2.4,  4.8,  6.0 
MHz,  although  the  4.8  MHz  is  used  for  normal  routine  operation.  Initially 
there  was  only  a  single  receiver  at  Westwood,  but  now  four  receiving  stations 
(Westwood,  Catskill,  Thornhurst  and  Lebanon)  are  in  full  operation.  With  the 
equipment  now  available,  multiple-frequency  measurements  at  a  single  station 
are  workable  as  shown  in  a  later  section. 

Microbarograph  array 

The  microbarograph  array  at  present  in  operation  consists  of 
ten  long  period  microbarographs  in  a  roughly  rectangular  configuration.  The 
average  spacing  of  the  instruments  is  approximately  60  kilometers  and  the 
array  aperture  is  about  240  kilometers.  The  uata  is  transmitted  in  analog 
form  over  telephone  lines  and  recorded  digitally.  The  initial  operation  of 
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TABLE  2.1 

FIELD  STATION  COORDINATES 


Latitude  Longitude 

Station  N  W 


Sterling  Forest,  N.Y. 

410 

11* 

50” 

74° 

14’ 

40" 

Thomhurst,  Pa. 

410 

12’ 

30” 

75° 

36’ 

30" 

Catskill,  N.Y. 

42° 

10* 

59” 

73° 

55' 

56” 

Lebanon,  N.J. 

39° 

55* 

05” 

74° 

31’ 

18" 

Stony  Brook,  L.I. 

40 P 

54* 

30” 

73° 

07’ 

32" 

Lamont  (Palisades),  N.Y. 

41° 

00' 

19" 

73° 

54’ 

28" 

Nauraushaun,  N.Y. 

410 

03’ 

26" 

73° 

59’ 

35" 

Pawling,  N.Y. 

410 

33' 

45" 

73° 

36’ 

04” 

Monticello,  N.Y. 

410 

39* 

24" 

74° 

40’ 

55" 

Hackettstown,  N.J. 

40° 

51* 

07.5" 

74° 

49’ 

44” 

Wagner  College,  S.I. 

403 

36’ 

55" 

74° 

05’ 

46” 

Westwood,  N.J. 

40  3 

59’ 

00" 

74° 

02' 

00" 

Fig.  2.1.  Map  showing  the  geographic  distribution  of  the 
stations  of  the  Teledyne  Isotopes  array. 
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the  array,  after  its  transfer  to  Westwood,  was  delayed  for  several  months  due 
to  the  lengthy  AT&T  strike,  but  data  has  been  gathered  continuously  since 
January  of  this  year.  The  instruments  have  given  relatively  trouble  free 
performance,  and  apart  from  routine  maintenance  no  repairs  have  been  necessary. 
This  aspect  has  been  particularly  gratifying  since  it  encourages  the  belief 
that  the  microbarographs  could  continue  operating  with  minimum  maintenance  in 
the  coming  year. 

Magnetometer  Array 

Three  3-component  magnetometers,  located  at  Catskill,  Thorn- 
hurst  and  Lebanon,  have  also  been  in  operation.  The  equipment  has  proved 
very  reliable,  requiring  almost  no  maintenance.  No  development  of  any  kind 
has  been  made  on  this  equipment  in  the  past  year,  and  the  data  collected  has 
proved  very  useful  in  the  study  of  ionospheric  background  motions. 

2.1.2  Array  Logistics  and  Instrumentation 

Logistics 

The  three  major  field  stations  are  located  in  isolated  areas, 
free  of  electrical  interference,  on  state  owned  lands  in  New  York,  Pennsylvania 
and  New  Jersey.  Trailers  house  the  phase-path  doppler  sounders,  microbaro¬ 
graphs,  magnetometers  and  the  required  digitizing  electronics.  These  three 
field  stations  are  linked  by  telephone  line  to  the  master  control  station  at 
Westwood,  N.  J.,  where  the  incoming  data  is  processed  and  stored  on  incremental 
digital  magnetic  tape.  A  fourth  station  identical  to  the  ones  mentioned  above, 
except  for  the  telephone  link,  is  hard  wired  to  the  master  station  at  Westwood, 
N.  J. 
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In  addition  to  the  four  major  field  stations,  seven  auxiliary 
sensor  stations  consisting  of  a  single  mierobarograph  and  telemetering 
electronics  are  housed  in  municipal  buildings  such  as  fire  houses,  schools, 
etc.  These  seven  auxiliary  stations  transmit  their  data  in  FM  analog  form 
by  telephone  lines  to  the  nearest  field  station  for  digitization. 

Instrumentation 

Phase-path  doppler  sounder 

The  high  frequency  continuous  wave  doppler  array  consists 
of  a  transmitter  and  four  receiver  stations.  The  doppler  systems  have  been 
built  after  Davies  (1962)  with  some  modifications  to  meet  our  specifications. 

The  transmitter  is  simple  and  reliable.  It  is  a  Heathkit  DX-60B 
transmitter  modified  for  continuous  operation  and  increased  tuning  range. 
Frequency  stability  is  controlled  by  a  precise  crystal  oscillator  stable  to 

g 

1  part  in  10  per  24  hours.  A  single  dipole  antenna,  an  automatic  keyer  and 
a  telephone  line  control  link  complete  the  transmitter  system.  Figure  2.2 
shows  the  simplified  block  diagram  of  the  doppler  system. 

A  doppler  receiving  system  consists  of  two  receivers:  a  signal 
receiver  and  a  reference  receiver.  The  signal  receiver  is  a  Hammarlund  SP600 
which  receives  the  transmitted  signal  after  reflection  from  the  ionosphere. 
The  reference  receiver  is  an  in-house  design.  A  brief  description  of  the 
principle  of  operation  of  the  receiving  system  follows. 

The  received  ionospheric  signal  is  converted  to  an  intermediate 
frequency  which  becomes  the  I.F.  signal  (455  KHz  ±Af)  in  the  signal  receiver. 
The  reference  signal,  which  comes  from  a  precision  crystal  oscillator 
identical  to  the  one  used  at  the  transmitter,  is  converted  to  an  intermediate 
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frequency  (I.F.  reference  455  KHz)  in  the  BM7  mixer,  which  is  slaved  to  the 
HFO  (Local  Oscillator)  of  the  signal  receiver.  The  difference  between  the 
I.F.  signal  and  I.F.  reference  is  equal  to  the  doppler  frequency.  This 
difference  is  obtained  in  the  BM9  mixer.  In  practice  the  reference  frequency 
of  the  transmitter  and  of  the  reference  receiver  are  offset  by  2  or  3  Hz  to 
allow  for  a  decrease  in  the  doppler  frequency.  The  difference  frequency, 
which  is  in  the  order  of  a  few  Hz,  is  amplified  and  converted  into  a  DC  volt¬ 
age  in  the  vidar  323  frequency  to  voltage  converter.  The  output  of  the  f/V 
converter  is  then  passed  through  a  low  pass  RC  filter  with  a  time  constant 
T  =  RC,  prior  to  digitization.  That  is,  each  digitized  point  represents  the 
time  average 


Mt)=i  E  t)dt 

0  T  J.  i. 


(2.1) 


The  time  constant  t  has  been  selected,  according  to  the  sampling  rate,  to 
minimize  aliasing  effects.  During  routine  operations  a  sampling  rate  of  10 
seconds  is  used. 


The  general  specifications  of  the  phase-path  doppler  system  are 


given  below. 


Signal  Receiver  -  Hammarlund  SP  600 

Frequency  band  of  interest  1.35  to  7.4  MHz 


Selectivity 

I.F.  crystal  fitter 

Sensitivity 


200  Hz 
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Reference  Receiver  -  in  house  design 

Input  frequency  range  0.05  MHz  to  12  Mllz 

Output  frequency  range  1  KHz  to  1  MHz 

Maximum  input  level  1  V  RMS 

Stability  is  referenced  to  precision  crystal 
oscillator  1  part  in  lO^/day 

Mixer  -  Spectran  Balanced  Ring  Modulator  S5-BM9 
Input  frequency  range  5  KHz  to  500  KHz 

Output  frequency  range  0  to  495  KHz 

Input  levels  1  V  RMS 

System  Noise  Level  and  Stability 

The  system  noise  level  is  taken  to  be  0.1  Hz.  The 
results  of  a  test  to  measure  system  noise  are  shown  in  Figure  2.3.  To  per¬ 
form  this  test  the  input  was  shorted  and  the  output  was  recorded  for  several 
days.  It  can  be  seen  that  over  any  1-hour  interval  the  signal  deviates  from 

its  mean  by  an  amount  less  than  or  equal  to  0.02  Hz  and  that  the  signal  drift 

over  any  24-hour  interval  is  approximately  0.08  Hz.  The  stability  of  the 
system  is  given  by  the  stability  of  the  precision  oscillator  which  is  1  part 

o  7 

in  10°  per  24  hours  plus  the  internal  noisp  or  about  1  part  in  10  . 

K j crobarographs 

Initially  the  gain  and  bandwidth  of  each  microbarograph 
were  adjusted  to  meet  selected  values  for  a  given  input  pressure  change.  Then 
two  or  more  microbarographs  were  set  to  measure  simultaneously  an  atmospheric 
signal  at  the  same  location  for  a  24-hour  period  and  their  outputs  were 
recorded  in  different  channels  of  the  tape.  Cross-spectrum  analysis  was 
applied  to  the  digital  records  of  the  different  pairs  tested.  In  this  way 
the  coherence  and  phase  response  of  each  instrument  were  obtained. 
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The  microbarographs  were  developed  and  built  at  Hudson  Laboratories 
of  Columbia  University,  and  they  have  been  described  by  Clay  and  Kraft  (1967)- 
Brief  specifications  are  listed  below. 


Frequency  response 

Sensitivity 

Range 


Flat  from  1  cph  to  180  cph 
3  pbar 
2500  pbar 


Magnetometers 

Two  kinds  of  magnetometers  are  used,  rubidium  vapor  Varian 

model  4938  and  Fluxgate  Honeywell  model  YK  10BI. 

Rubidium  magnetometer 
Frequency  response 
Sensitivity 


d.c.  to  20  Hz 
+0.01  gammas 


Fluxgate  magnetometer 
Frequency  response 
Sensitivity 


d.c.  to  10  Hz 
+0.1  gammas 


2.1.3  Array  Response  Calculations 

A  number  of  calculations  were  undertaken  to  compare  the  theoretical 
response  of  the  Isotopes  array  with  that  of  a  number  of  other  arrays;  a  10-element 
isometric  array,  a  15-element  isometric  array,  a  21-element  isometric  array  and 
the  LAMA  (Large  Aperture  Microbarograph  Array)  array.  The  isometric  arrays  have 
been  discussed  by  Haubrich  (1968)  and  would  appear  to  have  a  good  theoretical 
response.  The  spatial  arrangement  of  the  sensors  in  the  isometric  arrays  is 
given  by  Haubrich  and  the  aperture  was  chosen  to  be  225  km  so  that  it  would  be 
comparable  with  the  Isotopes  and  LAMA  arrays. 
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The  response  function  is,  of  course,  just  the  Fourier  transform  of  the 
array,  which  is  represented  by  a  set  of  delta  functions,  one  at  each  instrument 
position.  Thus  if  r^,  n  =  1,2,...,  N  are  the  positions  of  the  instruments, 


then , 


K(k)  =  [[  E  E  6(r-r)  6(r-f* )  e',k‘r  e  ^‘r  if/ 3r 

llii  n  m 

J J  n-1  m~l 


(2.2) 


K(k) 


N  N 
=  E  E 
n-1  ml 


(2.3) 


is  the  response  function  of  the  array.  This  is  just  the  response  to  a  wave 

of  infinite  wavelength.  Since  one  is  interested  in  the  response  to  a  wave  of 

finite  length  incident  from  a  given  direction,  i.e.,  to  a  wave  of  wave  number 

k  ,  the  function  calculated  was 
o 


K  (k-k  ) 
o 


N  N 

E  E 

n-1  m-1 


i( k— k  )*(r  -r  ) 
o'  o'  n  m 


(2.4) 


In  all  the  calculation  the  response  function  was  normalized  byK(o)  which 
can  be  shown  to  be 

K  (o)  =  N^,  if  N  is  odd,  (2.5) 

K  (o)  -  N(N+l),  if  N  is  even. 

Therefore  the  function  plotted  in  the  figures  is 

H(k-ko)  =  K(k-ko)/K(„)  .  (2.6) 

Figures 2 .4  through 2 .8  show  the  spatial  arrangement  of  the  Isotopes, 

TAMA,  and  three  isometric  arrays.  For  purposes  of  illustration  the  array 
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response  function  has  been  calculated  for  an  incoming  wave  coming  from  225° 
(southwest).  There  is  nothing  special  about  this  direction.  The  results  are 
representative  of  other  directions  of  wave  incidence  as  well. 

In  all  the  calculations  reported  below 

1*1  =  l*J  * 

The  azimuth  is  direction  in  which  the  array  was  steered,  i.e.,  the  direction 
of  k  . 

Figures  2.9  through  2 .12  show  the  response  of  the  three  isometric 
arrays  for  wavelengths  of  60,  100,  200  and  300  kilometers.  These  show,  as 
could  be  expected,  that  the  main  lobe  widens  as  the  wavelength  is  increased; 
that  there  are  a  number  of  side  lobes,  the  largest  about  0.25  and  0.30,  and 
that  the  21-element  isometric  array  has  the  lowest  side  lobes  for  a  given 
wavelength. 

Figures  2 .13  through  2 .16  show  the  response  of  the  21-e  lement  isome¬ 
tric  array,  the  LAMA  array  (20  elements)  and  the  Isotopes  array  (10  elements) 
for  wavelengths  of  60,  100,  200,  and  300  kilometers.  These  show  that  all  of 
these  arrays  have  main  lobes  of  comparable  width  and  comparable  side  lobes. 

From  an  examination  of  the  results  of  these  calculations  one  can 
conclude  that,  at  least  for  wavelengths  between  60  and  300  km,  the  response  of 
the  Isotopes  array  (10  instruments)  is  very  similar  to  that  of  the  LAMA  array 
(20  instruments)  and  that  of  a  21-instrument  isometric  array.  It  therefore 
follows  that  not  much  improvement  in  array  response  could  be  expected  by  adding 
new  instruments  or  changing  the  spatial  arrangement  of  the  present  array.  It 
was  therefore  decided  to  leave  the  Isotopes  array  unchanged. 


I 
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2.2  Event-Related  Pnonuworia 

2.2.1  Atmospheric  Nuclear  Tests 

During  the  operation  of  the  Toledyno  Isotopes  array,  atmospheric 
effects  due  to  the  Chinese  thermonuclear  explosion  of  2l>  September  196b  were 
detected . 

a)  A  signal  with  periods  between  2  and  5  minutes  arrived  at  about 
1830  G  M.T.  Beam  form ing  analyses  gave  a  phase  velocity  of 

320  +_  20  m/sec  and  an  azimuth  of  arrival  of  10°,  that  is.  travell¬ 
ing  along  the  short  great  circle  path  from  the  Lop  Nor  test  site  to 
the  array. 

b)  At  about  0018  G.M.T.  ori  30  September  1969,  a  wave  train  having 
periods  between  12  and  18  minutes  was  observed.  The  direction 
of  arrival  resulting  from  the  beainformiug  analyses  was  about 
180°,  that  is  coming  from  the  south  by  the  Jong  great  circle 
path,  with  a  phase  velocity  of  S80  +  20  in/ sec. 

The  zero  time  computed  from  the  short  period  <2  -  5  min)  arrival  was  08:50 

G.M.G.  20  minutes.  The  group  velocity  for  the  long  period  (.12  -  18  min) 

computed  from  the  arrival  time  was  540  +  20  m/sec. 

The  signal  to  noise  ratio  is  low  in  both  ranges  as  lias  been  reported 
earlier  by  Tolstoy  (I9«>7,  1968).  Thus,  filtering  and  beam  form  ing  were  necess¬ 
ary  to  extract  the  signal  from  the  background-  filtering  of  the  time  series 
was  carried  out  by  Iran  stunning  them  into  the  frequency  domain  with  a  fast 
Fourier  transform;  the  coefficients  of  the  unwanted  frequencies  were  set  equal 
to  zero,  and  the  remaining  Fourier  terms  were  transformed  hack  t.o  the  time  domain. 

Prior  to  beainforining  the  data  were  passed  tlmnigb  a  3  -  5  minute  filter 
in  the  case  of  the  short  period  arrival  and  a  12  -  18  minute  filter  for  the  long 
period  arrival.  These  filters  gave  the  sharpest  beam  form  ing  results  (Figures 
2.17  and  2.18) 
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The  beamforming  for  the  short  period  show  a  large  side  lobe.  This 
side  lobe  appears  to  be  a  characteristic  of  the  array  as  suggested  by  Figure  2.19 
which  is  the  theoretical  array  response  function  for  a  wavelength  equal  to  72  km 
corresponding  to  a  period  of  4  minutes  and  a  phase  velocity  of  300  m/sec. 

The  longer  period  beamforming  shows  a  broad  response  but  there  are  no 
side  lobes.  The  theoretical  response  for  the  long  period  shows  the  same  charac¬ 
teristics, 

A  more  quantitative  picture  can  be  obtained  by  calculating  for  each 
direction  a  quantity  which  is  proportional  to  the  overage  power 


(2.7) 


where  p  is  the  sum  of  the  amplitudes  of  the  signals  delayed  according  to  an 
assumed  phase  velocity,  and  divided  by  the  number  of  elements  of  the  array. 
Figure  2. 20 shows  the  result  of  applying  this  technique  to  the  time  interval 
1830  -  1850  Z  for  an  assumed  phase  velocity  of  300  m/sec.  If  one  lets  the 
phase  velocity  vary  in  the  range  <  V  <  one  obtains  Figure  2.21  which 
is  the  perspective  view  of  a  three-dimensional  surface.  The  value  of  P  in 
this  figure  is  along  the  z  axis;  the  phase  velocity  along  the  y  axis,  and 
the  direction  of  arrival  along  the  x  axis. 

Application  of  this  method  to  the  time  intervals  1810  -  1830  G.M.T. 
(prior  to  the  arrival),  1830  -  1850  G.M.T.  (containing  the  signal);  and  1850  - 
1910  G.M.T.  (after  the  arrival)  permits  a  comparison  between  the  background  and 
the  signal  for  the  short  period  arrival  (Figure  2.22).  If  we  now  pick  up  the 
phase  velocity  for  which  P  shows  the  maximum,  and  the  procedure  is  applied  to 
overlapping  time  intervals,  a  time  history  of  the  average  power  as  it  crosses 
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FIG.  2.20  Normalized  result  of  the  square  and  integration  technique 
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the  array  is  obtained.  Figure  2.23  displays  the  perspective  view  of  the  surface 
(P(0,t)  for  a  phase  velocity  of  300  m/sec  and  a  bandpass  filter  of  3  -  5  minutes, 
that  is,  for  wavelengths  between  54  and  90  km  and  covering  the  time  interval 
1750  -  1910  G.M.T.  The  integration  window  was  20  minutes  long  and  the  overlapp¬ 
ing  was  18  minutes. 

Figure  2 .24  shows  the  result  of  applying  the  above  procedure  to  the 
long  period  arrival  for  a  phase  velocity  of  540  m/sec  and  a  bandpass  filter 
of  12  -  18  minutes.  The  integration  window  in  this  case  was  60  minutes  and  the 
overlapping  was  54  minutes. 

The  amplitude  of  the  pressure  fluctuation  of  the  long  period  wave 
was  52  iibars  (peak  to  trough).  The  short  period  wave  had  an  amplitude  of 
24  g,bars  (peak  to  trough). 

A  comparison  between  those  observations  and  the  observations  of 
Tolstoy  of  the  nuclear  events  of  Lop  Nor  on  17  June  1967  and  the  French  test 
in  the  South  Pacific  on  24  August  1968  (Table  2.2)  permits  us  to  draw  some 
conclusions! 

i)  Short  period  waves.  These  waves  arrive  with  group  velocities  of  the  order 
of  300  m/sec.  The  center  frequency  of  the  three  events  listed  in  the  table 
takes  different  values  in  each  case.  These  values  range  from  0.002  Hz  (8  minutes) 
to  0.0067  Hz  (2.5  minutes);  that  is,  the  bandwidth  of  these  arrivals  is 

0.00167  Hz  -  0.0167  Hz  or  ]  minute  -  10  minutes.  The  average  phase  velocities 
are  300  _+  10  m/sec  and  the  wavelengths  range  from  J8  km  to  D80  km.  The  arrival 
corresponding  to  the  long  great  circcl  path  cannot  be  recognized. 

ii)  Long  period  waves.  The  group  velocities  of  these  waves  are  of  the  order 

of  550  m/sec  arid  the  phase  velocities  of  the  order  of  570  30  m/sec.  The  cen¬ 
ter  frequencies  also  vary  from  event  to  event  as  in  the  short  period  ease.  For 
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FIG.  2.23  Result  of  square  and  integration  for  a  fixed  phase  velocity 
V  =  300  m.  sec-1  using  overlapping  time  windows.  The  data 
corresponds  to  the  short  period  arrival. 
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the  three  events  compared  the  center  frequencies  are:  3.64  CPM  (lb. 5  minutes), 
4.62  CPM  (13  minutes),  and  4.0  CPM  (15  minutes)  respectively.  The  bandwidths 
go  from  0.000925  Hz  to  0.00138  Hz  (18  minutes  -  12  minutes).  Normally  both 
the  short  and  the  long  great  circle  path  arrivals  can  be  identified.  However, 
the  short  great  circle  path  arrival  corresponding  to  the  29  September  1969 
explosion  could  not  be  extracted  from  the  background  while  the  long-path 
arrival  required  very  little  filtering  before  it  could  be  identified. 

iii)  It  is  interesting  to  know  that  the  energy  of  the  explosion  Can  be  propa¬ 
gated  in  such  a  wide  range  of  periods.  It  is  also  important  to  notice  that 
the  events  compared  may  be  telling  us  something  about  the  variation  of  the 
structure  due  to  seasonal  change.  In  other  words,  the  events  occurred  during 
late  spring,  mid-summer  and  early  fall  respectively. 

iv)  The  lack  of  short-path  arrival  for  the  event  of  29  September  1969  could 
probably  be  explained  by  the  interaction  of  a  wind  structure  with  the  travell¬ 
ing  wave  as  it  propagated  through  the  polar  regions  in  the  northern  hemi¬ 
sphere  . 

TABLE  2.2 


17  June 

1967 

3 

< 

CM 

just  1968 

29  Sept  1969 

Short 

Period 

Range 

Long 

Period 

Range 

Short 

Period 

Range 

Long 

Period 

Range 

Short 

Period 

Range 

Long 

Period 

Range 

Period 

T  (min) 

5-10 

15-18 

2-3 

12-14 

2-5 

12-18 

Phase  Velocity 
V  (m/sec) 

300+10 

600+20 

300 

550+20 

300+10 

540+20 

Group  Velocity 
U  (m/sec) 

300 

550 

300 

550 

300 

540 
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2.2.2  SATURN -APOLLO  IX  Associated  Signals 

During  February  a  6.03  MHz  doppler  transmitter  was  installed  at 
Sterling  Forest.  This  signal  is  reflected  at  a  higher  altitude  than  the  4.8 
MHz  signal.  Both  the  6.03  MHz  and  4.8  MHz  doppler  ionosounders  recorded  an 
interesting  phenomenon  at  the  time  of  the  APOLLO  IX  launch  on  March  3.  The 
raw  data  from  the  4.8  MHz  and  6.03  MHz  doppler  are  shown  on  Figure  2.25. 

The  vertical  scale  is  1.0  Hz/inch  and  the  horizontal  scale  is  10  min/inch. 

The  sample  rate  was  10  points/minute.  The  starting  time  for  the  data  shown 
in  Figure  2.25  is  16:13:30  GMT  on  3  March  1969.  Tt  would  appear  that  between 
17:00  and  17:10  there  is  definite  change  in  the  character  of  the  signal  on 
both  channels. 

The  raw  data  was  passed  through  a  bandpass  (0.2  to  2-0  cycies/min) 
digital  filter.  The  filtered  data  for  the  4.8  MHz  and  6.03  MHz  dopplers  are 
plotted  in  Figure  2.26.  The  change  in  the  character  of  the  signals  is  even 
more  apparent  on  these  plots.  This  change  takes  place  between  17:05  and  17:15. 
The  amplitude  of  the  oscillations  increases  and  the  average  period  decreases 
to  about  1  to  2  minutes.  This  change  lasts  for  about  40  minutes. 

Cape  Kennedy  is  approximately  1500  kilometers  from  Westwood.  For 
a  signal  travelling  at  a  speed  of  500  m/sec,  the  travel  time  is  about  50 
minutes.  It  is  therefore  possible  that  this  signal  is  an  acoustic-gravity 
wave  which  was  generated  at. the  SATURN-AP0LL0  IX  launch  and  is  now  propagating 
in  the  atmosphere.  Since  the  planned  doppler  array  is  not  yet  functioning  it 
is  impossible  to  determine  the  direction  from  which  the  signal  came. 

Signals  with  the  same  characteristics  as  described  above  were 
recorded  during  the  SATURN-AP0LL0  XII  launch  on  14  November  1969.  At  that 
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FIG.  2.26  Filtered  doppler  records  corresponding  to  Saturn  Apollo  IX 
generated  signals.  Band  pass  filter  was  0.5  -  5  rain. 
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time  the  phase-path  sounder  array  consisted  of  four  elements  at  4.8  MHz  and 
one  element  at  6.03  MHz.  The  direction  of  arrival  obtained  from  the  array 
analyses  was  about  170°  measured  GW  from  geographic  north.  This  confirms  our 
early  hypothesis  about  the  origin  of  the  signals  on  3  March  1969. 

2.2.3  Correlation  of  Atmospheric  Pressure  Waves  with  Ionospheric 
Phase-path  Sounder  Signals _ 

During  the  summer  of  1968,  Hudson  Laboratories  of  Columbia  Univer¬ 
sity  was  operating,  north  of  New  York  City,  a  large  aperture  (several  hundred 
kilometer)  microbarograph  array  for  the  study  of  long-period  atmospheric 
pressure  waves.  The  array  (described  by  Herron  et  al.,  1969)  was  recording 
data  throughout  the  interval  during  which  there  were  five  high  energy 
disturbances  in  the  South  Pacific.  In  the  latter  part  of  the  summer,  a 
vertical  incidence  ionospheric  doppler  sounder  was  installed  at  the  center 
of  the  array  and  was  recording  during  the  event  of  24  August  1968.  The 
doppler  sounder  measured  the  rate  of  change  of  vertical  motion  of  electrons 
in  the  ionosphere.  A  cw  signal  (4.8  MHz)  was  transmitted  from  the  Sterling 
station  and  was  received  at  Hudson  Laboratories,  36  km  distant,  after  having 
been  reflected  from  the  lower  F  region  of  the  ionosphere  at  nearly  vertical 
incidence . 

After  the  project  was  transferred  to  Teledyne  Isotopes  the  analysis 
of  the  data  was  continued  and  the  principal  findings  are  described  below. 

Pressure  signals 

Tolstoy  and  Herron  (1970)  have  reported  the  existence  of  long- 
period  (10-25  min),  high-velocity  (600  m  sec  "*■)  atmospheric  pressure  signals 
identified  as  surface  gravity  waves.  Generated  by  high  energy  events  and 
detectable,  thus  far,  only  through  the  signal-to-noise  enhancement  of  array 
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processing  techniques,  the  waves  have  been  observed  travelling  the  short  and 
long  great-circle  paths  away  from  the  sites  of  the  events.  The  pressure  signal 
arrived  from  an  azimuth  of  230°  at  about  2300Z  on  24  August  1%8. 

Phase-path  sounder  signals 

At  the  time  at  which  the  pressure  signal  was  identified,  the 
phase-path  records  showed  a  group  of  sinusoidal  oscillations  with  periods 
between  10-15  minutes. 

Figure  2.27  shows  the  delayed  and  summed  pressure  trace  for 
230°  azimuth  above  the  filtered  phase-path  doppler  signal. 

Several  days  of  phase-path  sounder  records  before  and  after  the 
arrival  were  bandpass- filtered  (Figure  2.28).  It  must  be  noted  that  the  signal 
at  2300Z  on  24  August  1968  has  the  largest  amplitude. 

It  is  suggested  here  that  the  observed  phase-path  variation  is 
the  result  of  vertical  changes,  in  the  reflecting  region  of  the  ionosphere,  due 
to  the  passage  of  the  atmospheric  gravity  wave  generated  by  the  nuclear  event 
at  the  South  Pacific  test  site.  The  estimated  displacement  of  the  ionization 
was  of  the  order  of  5  km  peak  to  trough. 

2.3  Background  Fi  Id 

2.3.1  Ionospheric  Background  Motions  -  Relevance  to  Noise  Field  Phenomena 
This  section  is  directed  to  the  study  of  ionospheric  motions  with 
time  scales  of  the  order  of  minutes.  The  observations  have  been  made  with  phase- 
path  doppler  sounders  which  respond  to  vertical  motions  of  ionization. 

The  importance  of  this  study  rests  on  the  information  that  can  be 
obtained  about  the  background  of  ionospheric  motions  with  respect  to  the  driv¬ 
ing  mechanism  and  to  the  variations  of  the  structure.  From  the  theoretical 
point  of  view  it  is  interesting  because  it  provides  the  orders  of  magnitude  of 


2-42 


FIG.  2.27  Beamformed  pressure  trace  for  230°  azimuth  is  plotted  above 
filtered  doppler  record.  The  data  corresponds  to  the  French 
test  on  August  24,  1968. 
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the  motions  characterizing  the  background.  The  variation  of  the  structure  is 
an  important  clue  to  the  arrival  of  atmospheric  waves  originating  in  nuclear 
detonations.  From  the  point  of  view  of  communications  this  knowledge  is 
necessary  in  applications  such  as  over-the-horizon  radars  which  depend  on 
reflections  from  the  ionosphere  for  the  detection  of  missiles  and  aircraft. 

Although  several  mechanisms  may  contribute  to  the  motion  of  electrons 
in  the  ionosphere,  two  broad  categories  can  be  distinguished,  a)  Those  caused 
by  electrodynamic  forces  such  as  hydromagnetic  waves,  plasma  instabilities, 
and  electric  fields,  b)  Those  caused  by  motion  of  the  neutral  (un- ionized) 
component  of  the  atmosphere  (winds,  turbulence,  gravity  waves),  and  coupled  to 
electron  density  perturbations. 

Two  types  of  observations  are  discussed  in  this  section.  One  type 
requires  a  hydromagnetic  wave  interpretation  and  thus  belongs  to  the  first 
category.  The  other  type  can  be  explained  as  interaction  of  atmospheric 
gravity  waves  with  ionization  and  therefore  it  belongs  to  the  second  category. 

Mid-latitude  ionospheric  characteristics 

The  ionosphere  extends  upwards  from  60  to  500  km.  The  levels 
below  90  km  constitute  the  ionospheric  D  region.  Between  90  km  and  130  km  is 
the  E  region  where  a  well-defined  layer  of  ionization  is  formed  during  normal 
daytime  conditions.  Above  130  km  two  layers,  FI  and  F2,  are  formed  during 
daytime  hours;  this  region  is  called  the  F  region.  The  ionization  density 
normally  reaches  a  peak  in  the  F  region,  at  heights  of  250  to  5oo  km.  These 
heights  are  of  course  average  heights  since  they  vary  with  solar  time  and 
geomagnetic  location.  The  ionization  of  the  D,  E  and  F  layers  is  produced 
mainly  by  the  UV  radiation  from  the  sun  and  disappears  at  sunset  due  to  rapid 
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recombination  rates.  The  F2  layer  is  produced  partly  by  photoionization  and 
partly  by  collision  of  neutral  particles  with  energetic  particles.  The 
recombination  rates  are  slow  and  consequently  this  layer  is  present  through 
the  night. 

The  composition  of  the  ionosphere  varies  from  layer  to  layer.  Thus 
the  D  region  is  characterized  by  dominant  ions  of  N.^  and  0,^.  At;  E  region 
heights  0.^  ions  appear  to  be  dominant;  although  the  abundance  of  N.^  is  about 

4-  4- 

the  same  as  0,^,  the  high  recombination  rate  of  N,^  removes  the  N.^  ions  as 

rapidly  as  they  are  produced.  The  dominant  ion  at  F  region  levels  below  600 

4-  4  4- 

km  is  0  ,  although  a  very  small  fraction  of  the  ions  present  are  0.^  and 

ions . 

In  the  D  and  E  regions  collisions  of  electrons  with  neutral  particles 
predominate.  This  is  not  necessarily  the  case  in  the  F  region,  where  for  most 
purposes,  the  important  collisions  are  between  ions  and  neutral  particles. 
Mid-latitude  geomagnetic  fluctuations 

Rapid  magnetic  variations,  with  time  scales  varying  from  a  few 
tenths  of  a  second  to  many  hours,  and  amplitude  scales  varying  from  a  few 
hundredths  to  a  few  tens  of  gammas,  are  observed.  The  larger  amplitudes 
correspond  to  the  longer  period  fluctuations.  Geomagnetic  micropulsations 
are  classified  according  to  their  character  and  frequency  content  (Jacobs  et 
al.,  1664).  From  experimental  observations  two  main  classes  are  recognized: 
those  with  a  regular  and  continuous  character,  and  those  with  an  irregular  and 
impulsive  pattern.  They  are  divided  in  sub-groups  according  to  their  periods. 
The  regular  class  was  divided  into  five  period  ranges  and  designated  Pci 
through  Pc5.  The  irregular  class  was  divided  into  two  period  ranges  and  given 
the  designation  Pil  and  Pi2. 
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Notation 

Period  Range 

?  se 

Pci 

0.2 

to 

5 

Pc  2 

5 

to 

10 

Pc  3 

10 

to 

45 

Pc4 

45 

to 

150 

Pc5 

150 

to 

600 

Pil 

1 

to 

40 

Pi2 

40 

to 

150 

It  must  be  pointed  out  that  this  classification  is  a  guide  rather 
than  a  clear-cut  distinction. 

Geomagnetic  micropulsations  are  measured  with  a  variety  of  instru¬ 
ments:  nuclear  resonance  magnetometers,  fluxgate  magnetometers,  induction 
coils,  earth-current  meters. 

Each  type  of  micropulsation  shows  different  diurnal  and  seasonal 
variations,  thus  Pc  activity  shows  a  maximui  near  noontime  and  a  minimum 
during  winter  months.  Pi  activity  has  a  maximum  at  midnight  hours  but  does 
not  show  a  definite  seasonal  variation,  although  some  suggestions  of  an 
equinox  maxima  have  been  made. 

Previous  observations  of  ionospheric  motions 
In  order  to  put  the  studies  carried  out  under  this  contract  in 
proper  perspective,  it  is  necessary  to  review  briefly  the  previous  observa¬ 
tions  of  ionospheric  disturbances. 

As  early  as  1950,  Munro  in  Australia  considered  that  travelling 
disturbances  in  the  ionosphere  were  associated  with  disturbances  in  the 
neutral  atmosphere.  Martyn  (1950)  interpreted  the  wave-like  motions  observed 
by  Munro  (1950)  as  produced  by  internal  atmospheric  gravity  waves  that  provided 
continuously  rotating  cells,  of  neutral  gas,  in  the  F  region:  the  cells  were 
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moved  horizontally  by  a  wind  motion.  This  idea  was  however  abandoned  because 
of  the  direction  of  motion  required  by  the  model.  In  the  following  years 
other  mechanisms  involving  the  effects  of  dominant  electric  fields  or  hydro- 
magnetic  waves  were  proposed  (Martyn,  1955;  Dungey,  1955). 

By  the  end  of  the  1950's  several  causes  of  ionospheric  irregulari¬ 
ties  were  suggested.  Some  of  the  principal  ones  were:  variability  in 
ionization  rates,  mixing  by  air  currents,  recombination  and  diffusion  by 
winds,  hydromagnetic  waves,  turbulence,  gravity  waves.  Although  it  was 
recognized  that  fluid  motions  in  the  atmosphere  played  an  important  role,  the 
detailed  nature  of  this  role  was  a  matter  of  considerable  controversy. 

In  1957  at  the  URSI  General  Assembly  it  was  decided  that  a  joint 
study  by  a  group  of  ionospheric  physicists  and  fluid  dynamists  would  prove 
worthwhile.  A  symposium  was  then  organized  in  which  several  driving  mechanisms 
were  discussed.  Some  of  the  conclusions  of  the  symposium  were  (Bolgiano, 

1959): 

1)  Turbulence  is  a  common  occurrence  in  the  ionosphere,  at  least  to  an 
altitude  of  100  km.  This  was  established  from  the  rate  of  diffusion 
of  visible  meteor  trails,  which  is  several  orders  of  magnitude  higher 
than  that  attributable  to  molecular  diffusion  alone. 

2)  The  large-scale  anisotropic  motions  are  not  adequately  described  by 
the  theory  of  turbulence.  The  bulk  of  the  energy  probably  resides  in 
motions  more  properly  classed  as  semi-coherent  winds:  random  gravity 
waves,  thermal  winds.  The  intensity  of  fully  developed  turbulence 
may  be  no  more  than  1  percent  of  the  total  energy. 

3)  The  earth's  magnetic  field  has  negligible  direct  effect  on  turbulence 
but  may  affect  it  indirectly  via  the  dynamo  action. 

4)  Hydromagnetic  effects  probably  account  for  a  number  of  ionospheric 
phenomena . 
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In  1060  Hines  published  a  paper  in  which  he  explained  the  wind 
structures  revealed  by  the  meteor  trail  observations  of  Liller  and  Whipple 
(1954)  and  Greenhow  and  Neufeld  (1959)  in  terms  of  internal  atmospheric  gravity 
waves;  he  further  suggested  that  observations  of  drift  of  radio  diffraction 
patterns  produced  bv  moving  irregularities,  and  travelling  ionospheric  disturb¬ 
ances  could  also  be  explained  by  these  waves.  Ilieminger  and  Kohl  (1  962) 
suggested  that  the  500  m/sec  velocity  found  in  their  observations  were  sound 
waves  travelling  in  the  ionosphere,  where  the  speed  of  sound  is  higher.  Hines 
(1964)  pointed  out  that  internal  gravity  waves  can  account  for  a  wide  variety 
of  mid-latitude  irregularities,  suggesting  that  the  observed  propagation  to 
long  distances  without  attenuation  was  probably  due  to  ducting.  Wickersham 
(1965)  interpreted  the  natural  travelling  ionospheric  disturbances,  observed 
by  Munro  (1958)  and  others, as  ducted  acoustic-gravity  waves.  The  ionospheric 
effects  associated  with  the  Alaskan  earthquake  in  1964,  observed  by  Davies  and 
Baker  (1965),  were  explained  by  Row  (1967)  as  gravity  waves.  Georges  (1967) 
examined  single-frequency  vertical  incidence  CW  phase-path  doppler  data  from 
Boulder  with  regard  to  spectral  characteristics ;  he  suggests  that  the  longer- 
period  frequency  fluctuations  may  be  due  to  the  interaction  of  superimpose! 
gravity  waves  with  the  ionosphere  and  that  the  diurnal  variation  of  the 
dominant  period  is  an  observational  selection  due  to  the  change  in  radio 
reflection  height.  Titheridge  (1968)  interpreted  his  observations  in  terms 
of  gravity  waves,  pointing  out  that  the  periods  of  the  disturbances  were  near 
the  period  of  thermobaric  resonance.  Klostermeyer  (1969)  suggested,  from 
observations  of  atmospheric  gravity  waves  in  ionosonde  records,  that  the 
source  of  the  waves  was  in  the  troposphere.  Herron  and  Montes  (1970)  have 
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explained  their  observation  of  a  doppler  signal  connected  with  the  high  energy 
event  of  24  August  1968  as  a  surface  gravity  wave  mode.  Thome  and  Rao  (1969) 
concluded  from  their  observations  of  TIDs  over  Puerto  Rico,  that  the  principal 
effect  of  the  gravity  wave  motion  of  the  neutral  gas  is  to  impart  its  motion, 
parallel  to  the  earth's  magnetic  field,  to  the  ions  through  collisional  inter¬ 
action,  and  that  as  a  result  of  the  interplay  of  a  number  of  factors,  certain 
Fourier  components  of  the  permitted  gravity  wave  spectrum  are  more  effective 
than  others  in  producing  the  disturbances. 

Previous  observations  of  phase-path  variations  associated  with 

geomagnetic  activity _ 

Rishbeth  and  Garriot  (1964)  discussed  two  possible  mechanisms 
for  the  association  of  phase-path  fluctuations  and  geomagnetic  micropulsations 
found  by  Chan  et  al.  (1962).  The  first  mechanism  consists  of  a  hydromagnetic 
wave  travelling  along  the  geomagnetic  field  lines,  which  are  distorted  into 
a  sinusoidal  form;  the  wavelengths  of  the  displacements  are  greater  than  the 
vertical  extent  of  the  ionosphere.  The  second  mechanism  considers  that  the 
micropulsations  are  due  to  alternating  electromotive  forces  generated  by  the 
dynamo  action  of  the  oscillations  of  the  neutral  gas.  Unfortunately  the  iono¬ 
spheric  data  of  Chan  was  taken  at  oblique  incidence  including  several  hops  of 
che  HF  propagation  paths;  also  the  distance  between  reflection  points  and 
magnetometers  was  greater  than  1000  km  and  they  could  not  distinguish  between 
the  two  effects.  Jacobs  and  Watanabe  (1966)  made  a  more  detailed  treatment  of 
the  alternating  current  model.  They  concluded  that,  on  those  occasions  where 
doppler  shifts  of  about  1  Hz  corresponding  to  geomagnetic  micropulsations  with 
periods  of  50  to  90  sec  and  magnitudes  of  5  to  10  gammas,  the  micropulsations 
could  not  be  generated  by  alternating  ionospheric  currents.  Boyd  and  Duffus 
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(1969)  designed  an  extensive  experiment  to  observe  ionospheric  motions  asso¬ 
ciated  with  geomagnetic  micropulsations;  unfortunately  they  failed  to  do  any 
cross-spectral  analyses  of  the  records,  thus  phase  and  coherence  information 
were  not  obtained  by  these  workers,  and  their  conclusion  was  that  the  two 
most  likely  mechanisms  for  the  association  of  micropulsations  and  phase-path 
fluctuations  were:  (i)  hydromagnetic  waves  and  (ii)  neutral-gas  waves 
modulating  ionospheric  electrojots. 

Simultaneous  observations  of  phase-path  and  magnetic  variations 

During  the  month  of  August  1968,  as  part  of  the  atmospheric  wave 
propagation  program  of  Hudson  Laboratories  of  Columbia  University,  a  three- 
station  3-component  array  of  magnetometers  was  operated  in  conjunction  with  a 
vertical  incidence  phase-path  doppler  sounder,  and  an  array  of  twelve  micro¬ 
barographs  described  by  Herron  et  al-  (1969).  The  distribution  of  the  stations 
is  shown  in  Figure  2.1.  The  operation  of  the  array  continued  during  the  month 
of  September  after  the  project  was  transferred  to  Teledyne  Isotopes. 

The  frequency  used  in  the  doppler  sounder  was  4.824  MHz.  The  average 
reflection  height  for  this  frequency  is  of  the  order  of  150  km  during  daytime 
hours  and  200  km  during  nighttime  hours.  The  transmitter  was  located  at  Ster¬ 
ling  Forest,  and  the  receiver  at  Dobbs  Ferry,  New  York.  The  separation  between 
these  two  locations  was  about  36  km;  for  a  reflection  height  of  200  km  the 
angle  of  incidence  would  be  in  the  order  of  5  degrees.  For  all  practical 
purposes  we  could  assume  that  the  operation  was  at  vertical  incidence.  The 
ground  wave  was  conveniently  eliminated  by  the  orientation  of  the  antennas  and 
the  shielding  effect  of  the  Ramapo  Mountains  which  are  located  between  the 
transmitting  and  receiving  sites.  The  sampling  rate  was  10  seconds;  the  analog 
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signal  was  passed  through  a  10-second  low  pass  RC  filter  prior  to  digitization 
to  avoid  aliasing  effects. 

Boyd  and  Duffus  (1969)  reported  previous  observations  of  this  sort. 
However  there  were  some  unanswered  questions  related  to  the  mechanism  causing 
the  coupling  between  the  magnetic  micropulsations  and  the  phase-path  changes. 
Having  microbarographic  data  available  as  well  as  an  array  of  magnetometers, 
we  thought  that  these  questions  could  be  answered. 

Data  analysis 

During  the  period  of  recording  we  found  16  events  (see  Table  2.3) 
in  which  an  almost  peak-to-peak  correlation  between  magnetometer  and  doppler 
data  exists.  Note  that  the  times  of  occurrence  are  before  or  near  local  mid¬ 
night.  The  results  of  the  analyses  of  a  typical  signal  are  described  in  this 
report.  Typical  amplitudes  are  1  Hz  for  the  doppler  signal  and  1  to  2  gammas 
for  the  magnetic  signals. 

Power  and  cross-spectrum  analyses  using  the  method  of  Blackman  and 
Tuckey  (1959)  were  applied  to  the  doppler  and  magnetic  signals.  The  estimates 
of  the  power  and  cross  spectral  densities  are  given  by  the  Fourier  transform 
of  the  auto  and  cross-correlation  functions  respectively. 
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TABLE  2.3 
Events  Recorded 


Date 
8/27/68 
8/27/68 
8/31/68 
9/  2/68 
9/  3/68 
9/  4/68 
9/  5/68 
9/17/68 
9/17/68 
9/18/68 
9/18/68 
9/19/68 
9/21/68 
9/26/68 
9/27/68 
9/29/68 


Universal  Time 

0300  - 

0600 

0930  - 

1000 

0315  - 

0415 

0557  - 

0657 

0143  - 

0243 

0245  - 

0345 

0130  - 

0215 

0105  - 

0230 

0400  - 

0422 

0415  - 

0515 

2348  - 

0010 

0345  - 

0355 

0420  - 

0540 

0835  - 

1035 

0350  - 

0580 

0100  - 

0300 
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T/2 

Pu(t)  =  ^  ^  J  fx(t)  f^t+^dt  (2.10) 

-T/2 

T/2 

P12(t)  =  T-"-  1  J  fl(t)  f2(t+T)dt  (2,11) 

-T/2 


and  f^(t)  and  fg(t)  are  the  two  time  series  under  analyses. 

Prior  to  the  power  and  cross-spectral  analyses  the  data  was  passed 
through  a  pre-whitening  filter  which  consisted  of  a  simple  derivative  filter. 
The  smoothed  spectral  estimates  were  obtained  using  a  "hanning  window"  and 
10%  lags.  All  filters  used  in  the  process,  including  the  electric  filters  of 
the  data  acquisition  system,  were  restored. 


The  coherency  estimates  were  found  from 


|  Y(itu)  | 


(2.12) 


where  |y(iou)  |  =  magnitude  of  the  complex  coherency. 

Another  feature  worthy  of  study  is  the  phase  relationship.  The 
cross-spectral  analyses  give  the  phase  information  in  the  form 

ae(«)  =  02(<d)  -  9lU)  -  tan'1  {  j-  (2.13) 

where 


e-^ou) 

=  phase  of  time  series  f-^(t) 

(2.14' 

02(uj) 

=  phase  of  time  series  f2(t) 

(2.35) 

In  the  cross-sr ectrum  figures  the  names  identifying  the  two  time 
series  under  analysis  are  printed  at  the  top  of  the  figure.  The  starting  and 
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ending  times  are  written  next.  The  coherence  is  plotted  in  a  linear  scale 
from  £  to  1.  The  phase  is  indicated  by  two  square  symbols  joined  by  a  verti¬ 
cal  line;  the  length  of  the  line  is  the  90%  confidence  limit  for  the  phase 
angles.  When  no  vertical  line  is  drawn  and  only  one  square  symbol  appears 
the  confidence  limit  has  exceeded  180  degrees  and  the  symbol  is  drawn  at  the 
value  of  the  phase  given  by  the  analysis.  The  curves  in  the  lower  part  of  the 
figure  represent  the  power  spectrum  of  each  time  series.  The  power  spectrum 
levels  have  been  multiplied  by  the  center  frequency  of  the  band  in  which  the 
spectral  estimate  is  calculated.  The  continuous  curve  corresponds  to  the  time 
series  whose  identification  name  is  listed  first,  and  the  curve  with  the 
crosses  corresponds  to  the  second  time  series.  Each  cross  represents  a  fre¬ 
quency  band . 

A  typical  event  is  shown. 

Event  of  21  September  1968 

The  records  for  the  time  interval  0320Z-0550Z,  corresponding  to 
the  three  magnetometers  at  Catskill  and  the  phase-path  doppler  sounder,  are 
shown  in  Figure  2.29.  The  records  for  Thornhurst  and  Lebanon  are  practically 
identical. 

A  group  of  pulsations  starting  at  about  0420Z  and  lasting  until  0540Z 
can  be  correlated  in  all  the  traces;  this  is  better  shown  in  Figure  2.30  which 
represents  the  same  traces  filtered  with  a  Fourier  band-pass  filter  (0.5  min  - 
3.5  min). 

The  power  spectral  results  for  three  consecutive  windows  of  three 
hours  each  indicated  that  at  the  time  of  the  event,  there  is  an  increase  in  the 
energy  density  of  the  doppler  fluctuations  and  of  the  geomagnetic  micropulsa¬ 
tions.  This  increase  in  the  energy  density  is  about  an  order  of  magnitude  and 
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FIG.  2.29  Unfiltered  traces  corresponding  to  the  magnetometers  at 
Catskill  and  the  phase-path  sounder  for  the  September  21 
1968  event. 


FIG.  2.30  Filtered  traces  corresponding  to  the  magnetometers  at  Catskill 
and  the  phase-path  sounder  for  the  September  21,  1968  event. 
The  filter  band  pass  was  0.5  -  3.5  minutes. 
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it  is  mainly  confined  in  the  frequency  band  0.008  Hz  -  0.018  Hz  and  it  drops 
after  the  event,  Figure  2.31.  Each  trace  has  been  shifted  one  unit  in  the 
vertical  scale  to  avoid  superposition. 

Gross-spectral  analyses  have  been  applied  to  the  pairs:  NS 
Catskill  -  doppler;  EW  Gatskill  -  doppler;  EW  Catskill  -  NS  Catskill;  Rb 
Gatskill  -  doppler;  NS  Catskill  -  Rb  Gatskill;  EW  Catskill  -  EW  Thornhurst; 

NS  Catskill  -  NS  Thornhurst;  EW  Gatskill  -  Rb  Catskill;  EW  Gatskill  -  EW 
Lebanon;  EW  Lebanon  -  EW  Thornhurst.  The  results  are  shown  in  Figures  2.32 
through  2.41,  and  they  will  be  discussed  in  terms  of  power,  coherence,  and 
phase  separately. 

Coherence .  The  coherence  between  magnetic  components,  of  the  same 
polarization,  at  different  stations  is  in  general  high  throughout  the  fre¬ 
quency  range  covered  by  the  cross-spectral  analyses  (Figures  2.37  and  2.40). 

The  pair  NS  Gatskill  -  doppler  showed  high  coherence  (>0.8)  in  the 
frequency  band  corresponding  to  periods  between  1  and  2  minutes  (Figure  2.32) 
The  coherence  for  the  EW  Gatskill  -  doppler  pair  was  rather  low,  except  at 
the  frequency  corresponding  to  1  minute  period  where  it  reaches  a  value  of 
about  0.6  (Figure  2 .33) .  The  Rb  Catskill  -  doppler  pair  followed  closely  the 
results  for  the  NS  Gatskill  -  doppler  pair  (Figure2.35) .  The  same  relation¬ 
ship  holds  between  the  magnetic  components  at  the  other  two  stations  and  the 
doppler.  The  EW-NS  pairs  at  each  station  showed  coherence  values  near  0.6 

at  periods  of  about  1  minute  (Figure  2.34);  this  is  the  same  type  of  relation 

that  exists  between  the  EW  -  doppler  pairs.  The  coherence  for  the  NS  -  Rb 

pair  showed  a  peak  between  1  and  2  minute  period  (Figure  2.36). 

Phase .  The  phase  relationships  are  as  follows.  The  NS  magnetic 
components  lag  the  doppler  signal  in  phase  by  about  140  degrees  at  a  period 
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FIG.  2.32 


Cross- spectrum  between  NS  Catskill  and  phase-path  doppler. 
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2.35  Cross-spectrum  between  Rb  Catskill  and  phase-path  doppler. 
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FIG.  2.36  Cross- spectrum  between  NS  Catskill  and  Rb  Catskill. 
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of  one  minute  decreasing  more  or  less  uniformly  to  about  116  degrees  at  2 
minute  period  (Figure  2.32).  Approximately  the  same  relationship  holds  for 
the  pair  Rb  -  doppler  (Figure  2.35).  The  phase  between  EW  Gatskill  and 
doppler  shows  considerable  scatter  for  periods  at  which  the  coherence  is  low, 
as  one  should  expect  in  theoretical  grounds,  since  coherence  is  a  measure  of 
the  phase  stability.  However,  for  periods  of  about  1  minute  the  phase  is 
relatively  stable  and  its  value  is  about  -15  degrees,  that  is  the  EW  compo¬ 
nents  lead  the  doppler  signal  in  phase  by  about  15  degrees  at  the  1  minute 
period  band  (Figure  2.33).  Important  phase  relationships  are  those  between 
magnetic  components  of  different  polarization  at  the  same  station;  and 
between  magnetic  components  of  the  same  polarization  at  different  stations. 

The  first  type  of  relation  gives  information  about  the  orientation  of  the 
disturbance  vector,  in  the  plane  containing  the  two  components  under  analysis; 
and  the  second  about  time  shift  between  stations  and  consequently  can  lead  to 
the  determination  of  apparent  phase  velocities.  The  EW  components  lead  the 
NS  and  Rb  components  by  157  degrees  and  140  degrees  respectively  at  periods 
of  1  minute.  The  NS  and  Rb  showed  phase  shifts  of  -20  degrees  for  2  minute 
period  and  11.7  degrees  for  1  minute  period  indicating  that  in  the  first  case 
the  Rb  lags  the  NS  and  in  the  second  case  the  Rb  leads  the  NS  component.  The 
phase  velocities,  estimated  from  the  phase  differences  between  magnetic  compo¬ 
nents  of  the  same  polarization  at  different  stations,  are  of  the  order  of  100 
km/sec  ±  50  km/sec.  The  large  uncertainty  in  the  phase  velocity  is  due  to  the 
small  number  of  data  points  used  in  the  computation  of  the  cross-spectra 
(Goodman,  1957). 
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2.37  Cross- spectrum  between  EW  Cat  skill  and  EW  Thomhurst. 
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FIG.  2.38 


Cross- spectrum  between  NS  Catskill  and  NS  Thomhurst. 
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'FIG.  2.39  Cross-spectm.n  betwsi  EWCatskill  and  RB  Catskill. 
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Interpretation  of  Results.  The  impulsive  character,  the  frequency 
constant,  and  the  times  of  occurrence  of  the  events  indicate  that  the  micro¬ 
pulsations  associated  with  the  phase-path  fluctuations  belong  to  the  Pi2 
type  (Jacobs  et  al.,  1965;  Christoffel  and  iinford,  1966;  Saito  and  Matushita, 
1968). 

The  lack  of  infrasonic  activity  associated  with  the  events  and  the 
high  values  of  the  phase  velocities  indicate  the  hydromagnetic  character  of 
this  type  of  micropulsations  in  agreement  with  recent  theoretical  work^ 
(Rostoker,  1967;  Saito  and  Matushita,  1968). 

The  coherence  results  show  that  the  disturbance  is  made  up  by  the 
superposition  of  several  modes  of  hydromagnetic  waves.  For  example  the 
event  of  September  21  has  one  mode  centered  at  a  frequency  corresponding  to 
a  period  of  1  minute  which  involves  all  three  magnetic  components;  the  other 
modes  correspond  to  periods  greater  than  1  minute  and  affect  only  the  NS  and 
Rb  components.  On  the  other  hand  the  event  of  September  38  showed  modes 
centered  at  1.14  minutes  and  2.0  minutes  respectively,  in  which  the  three 
magnetic  components  are  affected;  and  one  mode  centered  at  1.6  minutes 
affecting  only  the  NS  and  Rb  components.  The  magnitude  of  the  coherence 
between  different  wave  components  is  representative  of  the  degree  of  polari¬ 
zation  of  the  wave  in  the  plane  containing  the  pair  of  components  (Bom  and 
Wolf,  1965);  thus  the  mode  centered  at  1  minute  is  50%  polarized  in  the  EW  -  NS 
plane  and  98%  polarized  in  the  NS  -  Rb  plane;  similarly  the  modes  centered  at 
1.14  minute  and  2.0  minutes  are  50%  polarized  in  the  EW  -  NS  plane  and  about 
70%  in  the  NS  -  Rb  plane.  Th  modes  not  affecting  the  EW  have  degrees  of 
polarization  of  about  80%. 
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The  amplitudes  of  the  Rb  fluctuations  in  both  events  are  smaller 
than  those  of  the  EW  and  NS  components.  This  fact  together  with  the  high 
coherence  between  the  NS  and  Rb  signals  suggest  that  the  Rb  disturbance  is 
produced  mostly  by  induction  of  the  NS  signal  on  the  ground.  If  this  inter¬ 
pretation  is  granted,  one  can  postulate  that  the  magnetic  perturbations  are 
due  to  transverse  hydromagnetic  modes,  of  waves  propagating  along  the  field 
lines 

A  simple  quantitative  interpretation  can  be  obtained  from  the 
equations  relating  the  fluid  velocity  and  the  magnetic  field  disturbance 
of  a  transverse  wave  in  an  incompressible  fluid  of  infinite  conductivity. 
Following  Alfve*n  and  Falthammar  (1963),  consider  the  geometry  of  Figure  2.41 
and  the  MKS  equations 


St 


dz 


and 


whe  re 


V  -  velocity  of  the  fluid 
x 

magnetic  vector  of  the  wave 

Earth's  magnetic  field 
H  magnetic  permeability 

p  plasma  density 


(2.16) 


(2.17) 
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The  solution  of  equation  (2 17)  has  the  form 
sin  m(t  —  z/Vg) 

where 

bQ  =  magnitude  of  the  magnetic  vector 

Vg  =  BQ(^p)  '  is  the  phase  velocity  of  the  wave. 

Substitution  of  equation  (2.18)  into  equation  (2.20)  gives 

Vx  =  -  (g,p)-1^2  bQ  sin  u)(t  -  z/Va) 
and 

Vx  =  •  bx  (w)'1/2  • 


(2.18) 


(2.1^ 


(2  .20) 


Rishbeth  and  Garriot  (1964)  have  pointed  out  that  at  F  region  heights 
the  ion  neutral  collision  frequency  is  large  enough  to  modify  the  phase  velocity 
of  the  wave  and  suggested  the  use  of  a  complex  dispersion  equation  which  modi¬ 
fies  equation  (2.20)  to  give 


V 

x 


(up v  in) 


-1/2  -j3tt/4 


(2.21) 


The  inclusion  of  the  collision  frequency  is  the  same  as  taking  into  account 
the  conductivity  of  the  medium.  Dungey  (1967)  suggested  that  the  density  p 
can  be  replaced  by  the  electron  density  N  multiplied  by  the  mass  of  the 
dominant  ion  M^.  In  the  F  region  of  the  ionosphere  the  dominant  ions  are 
0  ,  and  the  electron  density  can  be  taken  from  the  plasma  frequency  corre¬ 
sponding  to  the  sounding  frequency.  Typical  values  for  the  parameters  entering 
equation  (2.21)  are 
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4tt  x  10" 7  MKS 

1.33  x  10-26  MKS 

11  3 

2.85  x  10  electrons/m 


V  in  = 


and  equation  (2.21)  projected  on  the  vertical  direction  becomes 

Vh,  =  V  COS  I  =  4.0  x  1010  b  T"1/2  ^  COS  I 

t  x  x 

taking  T  =  60  sec  and  k  I  -  2  x  10  9  equation  (2.22)  gives 


Kl  «  3  m / sec 


The  velocity  as  observed  on  the  phase  path  sounder  data  is 


given  by 


(2.22) 


lVobs  I  "  If  Af  (2-23> 

using  f  =  5.0  MHz  and  Af  «  0.5  Hz  one  obtains  l^0^s I  w  15  m/sec. 

Equation  (2.22)  predicts  a  phase  shift  of  135  degrees  between  the  velocity 

and  the  magnetic  disturbance.  The  observed  phase  shift  between  the 
doppler  and  the  NS,  as  was  pointed  out  earlier,  is  of  about  140  degrees, 
at  a  period  of  1  minute,  decreasing  as  the  period  increases.  Alfven  and 
Falthammar  (1963)  have  showed  that  for  simple  transverse  waves  when  the  con¬ 
ductivity  is  taken  into  account  the  phase  shift  0  is  proportional  to  the 
inverse  of  the  period 

0  a  j  (2.24) 


I 
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Using  this  relation  one  could  explain  the  decrease  in  phase  between  the 
phase-path  fluctuations  and  the  NS  component  as  the  period  of  the  disturb¬ 
ance  increases. 

The  difference  between  predicted  and  observed  amplitudes  is  proba¬ 
bly  due  to  attenuation  of  the  magnetic  disturbance  when  traversing  the 
ionosphere  at  E  region  levels  where  the  conductivity  is  higher. 

Jacobs  and  Watanbe  (1966)  have  pointed  out  that  the  association 
between  phase-path  variations  and  geomagnetic  micropulsations  cannot  be  due 
to  ionospheric  currents  if  the  magnetic  disturbance  and  the  doppler  shift 
have  the  same  order  of  magnitude  and  periods  of  the  order  of  100  sec.  This 
theoretical  result  favors  our  interpretation  that  the  waves  are  of  hydro- 
magnetic  origin.  Boyd  and  Duffus  (1969)  suggested  two  general  mechanisms 
for  the  association:  i)  hydromagnetic  waves,  ii)  neutral  gas  waves; 
however  they  have  not  proposed  any  coupling  mechanism  which  can  explain  the 
amplitude,  phase  and  coherence  between  the  parameters  measured. 

We  can  conclude  then  that: 

1)  The  association  between  geopiagnetic  micropulsations  and  phase-path 
variations  of  the  type  investigated  in  this  work  can  be  explained  by  a 
superposition  of  various  inodes  of  transverse  hydromagnetic  waves,  at  least 
within  an  order  of  magnitude. 

2)  At  least  a  significant  portion  of  the  ionospheric  motions  in  this  part 
of  the  spectrum  have  a  hydromagnetic  origin  since  this  association  is  not  a 
rare  phenomenon. 
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Multiple- frequency  phase-path  doppler  observations 

Observations  using  three  sounding  frequencies  simultaneously,  at 
one  location,  were  carried  out  during  the  months  of  May  and  June,  1969. 

A  test  for  the  coherence  and  frequency  responses  of  the  three  systems  was 
run  by  recording  simultaneously  in  the  three  systems  the  same  ionospheric 
reflected  signal  at  the  same  location.  Cross  spectral  analyses  of  the 
records  indicated  that  for  all  practical  purposes  the  instruments  were 
identical. 

The  CW  radio  transmitters  were  located  at  Sterling  Forest,  N.  Y. 
and  the  signals  were  received  at  Westwood,  N.J.  The  distance  from  the 
transmitters  to  receivers  was  30  km,  and  for  practical  purposes  vertical 
incidence  could  be  assumed.  The  transmitting  and  receiving  antennas  con¬ 
sisted  of  simple  horizontally  polarized  half-wave  dipoles.  The  frequencies 
used  were  2.412  MHz,  4.824  MHz  and  6.030  MHz;  the  average  reflection  heights 
for  these  frequencies  were  of  the  order  of  130,  170  and  230  km  respectively. 

The  variations  of  the  received  frequencies  as  a  function  of  time  were  recorded 
in  digital  form  on  an  incremental  digital  magnetic  tape.  Since  the  periods  of 
interest  were  of  the  order  of  minutes,  a  sampling  rate  of  10  seconds  was  con¬ 
sidered  convenient.  As  before,  a  low  pass  electronic  filter  was  used  to  mini¬ 
mize  aliasing  effects.  Data  were  recorded  mainly  during  part  of  the  month  of 
June,  although  some  data  using  two  frequencies  were  collected  in  May.  In 
some  nights  the  critical  frequency  of  the  ionosphere,  that  is,  the  maximum 
plasma  frequency,  was  smaller  than  6  MHz  for  several  hours,  and  consequently 
there  were  no  reflections  from  6.03  MHz  during  these  time  intervals;  the  usable 
data  for  these  nights  consisted  only  of  two  frequencies. 
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Since  in  this  experiment  one  important  parameter  was  the  reflection 
height,  the  heights  for  the  three  sounding  frequencies  were  obtained  from  regu¬ 
lar  ionograms.  The  relationship  between  true  height  of  reflection  and  virtual 
height  as  obtained  from  ionograms  is  described  next. 

Relationship  between  virtual  height  and  true  height 

Ionospheric  heights  of  reflection  as  a  function  of  frequency  are 
routinely  obtained  at  a  national  network  using  ionosondes.  The  ionosonde 
data  is  in  the  form  of  ionograms  which  are  records  of  "virtual  height"  versus 
sounding  frequency.  The  virtual  height  h'  is  the  height  at  which  a  radio 
pulse  would  be  reflected  if  it  would  travel  at  the  speed  of  light  in  free 
space,  at  all  times  during  its  way  to  the  ionosphere  and  back. 

h'  =  -|  ct  (2.25) 


where  £  is  the  speed  of  light  in  free  space  and  t  the  time  it  takes  the 
pulse  to  go  to  the  ionosphere  and  return  after  reflection.  However  since 
the  velocity  of  the  pulse  is  the  group  velocity  U 


where  hr  is  the  true  reflection  height,  and  ^  is  the  group  refractive 
index.  Equation  (2  .26)  can  be  inverted  analytically  using  Abel's 
integral  equation,  when  certain  assumptions  such  as  the  absence  of  magnetic 
field  and  electron  collisions,  are  made  about  the  properties  of  the  media 
(Davies,  1969). 


(2.26) 


h(f  )  = 
'  v 


fv 

2  f  h'(f)df 
"  Jo  Vf2  -  f2 


(2.27) 
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For  this  work  the  magnetic  field  and  electron  collisions  were  taken  into 
account  and  a  numerical  true  height  analyses  (Titheridge,  1967)  was  used. 
The  need  for  the  true  height  analysis  becomes  apparent  if  one  looks  at 
Figure  2.42  which  shows  a  typical  ionogram  and  the  corresponding  true 
height  analysis. 

Ionograms  for  Hanover,  N.H.  were  made  available  to  us  by  Prof. 
M.G.  Morgan  of  the  Radio  Physics  Laboratory  of  Dartmouth  College.  The 
heights  of  reflection  for  each  power  spectrum  were  obtained  by  averaging 
the  heights  of  hourly  ionograms. 

Computation  of  N(z) 

The  distribution  of  the  Vaislila  frequency  N(z)  with  height  was 
obtained  from  Tolstoy  (1967).  The  temperature  distribution  was  computed 
using  the  empirical  formulas  from  the  U.S.  Standard  Atmosphere  Supplements, 
1966. 

T(z)  =  T-  -  (T»  -  T120)  *  exp  [-S(z  -  120)] 

where 

S  =  0.0291  •  exp  (— q2/2 ) 

and 

q  =  (T*  -  800)  [750  +  1.722  x  10-4  x  (T®  -  800)2]-1 

T12()  =  temperature  at  120  km 

=  temperature  at  the  top  of  the  ionosphere 
The  values  used  for  T-^q  and  T®  were  327°K  and  200°K  respectively  which  are 
the  suggested  values  for  Summer  time. 


(2.28) 


(2.29) 


(2.30) 
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FIG.  2.42  Typical  ionogram  showing  the  true  height  analysis 


Preliminary  Results 


Power  spectra  of  the  phase-path  sounder  records  showed  systematic 
peaks  and  slope  breaks,  which  appear  to  shift  towards  longer  periods  for 
greater  reflection  heights  (larger  sounding  frequencies).  Examples  of  this 
kind  of  behavior  are  illustrated  in  Figures 2. 43  and  2.44). 

The  distribution  of  the  period  corresponding  to  the  main  spectral 
peak  with  height  (sounding  frequency)  shows  a  pattern  that  resembles  the  dis¬ 
tribution  of  the  VSisalS  frequency  for  an  adiabatic  atmosphere  (Fig.  2.45). 
This  work  is  continuing  in  the  present  contract  and  new  results  will  be  re¬ 
ported  at  a  later  date. 

Conclusions 

The  main  conclusions  of  this  section  are: 

1)  The  geomagnetic  micropulsations  associated  with  phase-path  variations 
in  the  period  range  0.5  minutes  to  3  minutes  are  of  the  Pi2  type  and 
consist  of  the  superposition  of  various  modes  of  transverse  hydro- 
magnetic  waves.  The  estimated  phase  velocities  are  in  agreement  with 
values  of  the  Alfven  velocity  at  ionospheric  heights.  These  hydro- 
magnetic  waves  modulate  the  ionospheric  plasma  producing  the  observed 
phase-path  variations. 

2)  The  spectral  peaks  of  phase-path  doppler  records  occurring  at  periods 
greater  than  5  minutes  probably  correspond  to  internal  gravity  wave 
activity. 

3)  Inclusion  of  viscosity  into  the  theory  of  internal  gravity  wave  propa¬ 
gation  will  allow  a  quantitative  interpretation  of  the  phase-path 
variation  spectra,  with  the  possibility  of  deriving  the  parameters  of 
the  neutral  gas  structure  at  ionospheric  levels. 
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FIG.  2.44  Power  spectra  of  phase-p*th  doppler  records  for  6.030, 
4.824  and  2.412  MHz  corresponding  to  the  time  interval 
1810  -  2410  G.M.T.  on  June  10,  1969. 
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2.3.2  Microbarograph  signal  and  noise  fields 

Acoustic-gravity  waves  generated  in  the  atmosphere  by  nuclear  explo¬ 
sions  have  been  recorded  and  studied  for  many  years.  Many  observational 
characteristics  of  these  signals  have  been  reported  and  various  theoretical 
studies  of  source  and  propagation  have  been  conducted.  One  property  of  the 
signals,  however,  that  has'  not  been  quantitatively  studied  is  the  spatial 
coherence. 

This  section  presents  information  regarding  the  similarity  (i.e. 
coherence)  of  the  signals  and  of  atmospheric  pressure  background  noise  over 
distances  of  several  kilometers  to  about  200  kilometers.  Such  information  is 
required  to  design  an  array  of  microbarographs  for  the  detection  of  the  signals 
and  is  required  in  choosing  a  processing  technique  for  the  array  data.  The  spatial 
coherence  has  not  been  previously  studied  due  to  a  lack  of  adequate  data,  i.e. 
digitized  data  from  a  sufficiently  large  and  dense  array  of  microbarographs • 

The  data  is  now  available  from  the  Hudson  Laboratories'  large  aperture  micro¬ 
barograph  array.  The  array  was  recording  during  the  French  thermonuclear  test 
of  24  August  1968.  The  spatial  coherencies  of  acoustic-gravity  waves  from  that 
test  have  been  determined  by  cross-spectral  analyses  of  many  pairs  of  stations 
of  the  array,  at  various  separations.  For  comparison  with  the  signal  coheren¬ 
cies,  an  analysis  of  background  noise  coherencies  has  been  performed. 

Signal  coherence 

The  24  August  1968  French  nuclear  test  was  a  megaton  range  explosion 
that  generated  acoustic-gravity  waves  clearly  recorded  on  the  Hudson  Laboratories 
array.  Various  modes  of  propagation  were  recorded  in  the  1-  to  5-minute  period 
band.  The  separate  modes  were  found  by  passing  the  data  through  successive  1/2 
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octave  digital  filters.  The  fundamental  (3-5  minute  period)  and  the  first 
acoustic  (2-3  minute)  had  the  highest  signal  to  noise  ratios  and  were  chosen 
for  a  study  of  their  spatial  coherence.  Cross-spectra]  analyses  were  per fanned 
between  almost  all  combinations  of  pairs  of  stations  in  the  array  giving  coher¬ 
encies  for  many  station  separations. 

The  signals  from  the  French  test  site  in  the  South  Pacific  crossed 
the  array  travelling  into  50°  azimuth.  Coherencies  for  station  pairs  aligned 
+  25°  from  the  direction  of  signal  travel  (that  is,  50°  +_  25°)  are  shown  in 
Figure  2.46  for  the  first  acoustic  inode  (2-3  minute  period),  us  a  function  of 
station  separation.  Also  shown  are  coherencies  for  pairs  aligned  _+  25°  from 
the  normal  to  the  direction  of  travel.  It  is  clear  that  the  transverse  coher¬ 
encies  are  higher  than  the  parallel  coherencies  at  large  separations,  that  is, 
the  waves  are  long-crested.  All  coherencies  (signal  or  noised  of  course, 
approach  unity  at  small  separations. 

A  similar  analysis  was  performed  on  the  fundamental  acoustic  mode 
(3-5  minute  period).  Figure 2.47  shows  that  the  transverse  coherence  is  again 
greater  than  the  parallel  coherence.  A  study  of  the  spatial  coherence  of 
higher  acoustic  modes  and  of  the  gravity  modes  is  in  progress  and  will  be 
includ  I  in  the  final  copy  of  this  report. 

Noise  coherence 

The  data  for  noise  coherencies  were  chosen  from  two  intervals  of 
time.  One  interval  (l  hour  long)  was  chosen  two  hours  prior  to  the  signal 
arrival.  Tills  data  gave  the  noise  coherencies  (2-3  minute  band)  at  the 
larger  station  separations  of  Figure  2.48.  The  other  interval  of  time  was 
a  one-hour  sample  several  months  before  signal  arrival,  when  the  small  array 
in  the  insert  of  Figure  2.1  was  operating.  The  data  from  this  small  array 
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provided  the  coherency  figures  at  the  small  separations  (less  than  14  km)  of 
Figure  2.48.  We  see  that  over  a  wide  range  of  separations  the  coherencies 
are  almost  all  below  the  95%  confidence  level. 

In  order  to  find  the  short  period  limit  at  which  the  noise  becomes 
coherent,  24  cross-spectral  measurements  were  performed  between  stations  of 
the  small  array  and  mean  coherencies  were  computed  for  separations  of  1.4 
to  14  km.  Figure  2.49  shows  mean  coherencies  at  32,  4,  and  3-minute  periods. 
At  these  periods  the  coherencies  drop  to  a  level  of  "random"  coherence  in 
several  kilometers. 

Results 

Figure  2.50  summarizes  the  results  of  the  signal-  and  noise-coherenct 
study.  The  coherence  of  acoustic  signals  is  well  above  that  of  the  background 
noise  at  separations  greater  than  3  or  4  kilometers.  The  signal  coherence 
remains  quite  high  out  to  separations  of  50  km  (or  even  more  for  the  transverse 
coherencies) . 

It  is  apparent  that  even  for  arrays  of  many  tens  of  kilometers  aper¬ 
ture,  the  signal  coherence  is  high  enough  to  beamform  with  the  entire  array. 
Beamfortning  on  2-5  minute  periods  was,  in  fact,  quite  successful  over  the 
200  km-aperture  Hudson  Laboratories  array. 

Group  velocity  cutves  of  long  period  atmospheric  waves 

Lo'»g  period  (10-15  minutes),  high  velocity  (600  in/sec)  atmospheric 
pressure  waves  have  been  observed  on  the  Hudson  Laboratories  large  aperture 
microbarograpli  array.  Through  beamforming  techniques,  the  signals  have  been 
observed  travelling  both  the  short  and  the  long  great  circle  pathi  away  from 
the  sites  of  nuclear  detonations.  Group  velocity  curves  of  these  signals 
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have  been  determined  by  the  method  of  plotting  arrival  times  of  peaks  and 
troughs  of  the  wavetrain  versus  the  peak  or  trough  number.  Figure  2.51 
shows  the  dispersion  curves  from  three  nuclear  detonations.  The  curves 
labeled  D  are  obtained  from  the  direct,  short-path  arrivals.  The  curves 
labeled  A  are  from  the  antipodal,  long-path  arrivals. 

As  seen,  the  long  periods  travel  with  higher  velocity.  This  type 
of  dispersion  is  common  to  internal  gravity  waves,  surface  gravity  waves  and 
the  long  period  branch  of  the  fundamental  acoustic  mode. 

It  is  interesting  to  note  that  the  velocities  are  computed  as 
higher  for  the  antipodal  arrival.  The  group  velocities  are  computed  as  the 
total  distance  travelled  divided  by  the  total  travel  time.-  The  higher  anti¬ 
podal  velocity  suggests  that  the  travel  times  used  were  too  high,  that  is, 
that  these  very  long  wavelength  signals  were  not  initiated  at  the  instant  of 
detonation.  Rather,  it  took  many  tens  cf  minutes  to  start  the  waves  moving. 
Using  overly-long  travel  times  (from  the  moment  of  detonation)  has  the  great¬ 
est  depressing  effect  on  the  short  path  velocities.  If  the  travel  times  are 
shortened  to  bring  the  short  and  long  path  velocities  into  coincidence,  the 
group  velocities  lie  in  the  600-700  m/sec  range. 

Thus  these  results  not  only  tell  us  the  range  of  the  group  veloci¬ 
ties  and  the  sense  of  dispersion,  but  also  indicate  that  the  source  mechanism 
requires  many  tens  of  minutes  to  initiate  such  long  period,  long  wavelength 
waves . 
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FIG.  2. SI  Group  velocity  dispersion  curves  for  nuclear  detonations. 
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3.  THEORETICAL  ATMOSPHERIC  PHYSICS 

3 .1  Wave  Propagation 

It  is  known  (Herron  and  Tolstoy,  1969)  that  internal  gravity  waves 
are  associated  with  the  jet  stream.  It  has  been  hypothesized  that  these 
internal  gravity  waves  are  gene  ated  by  a  non-linear  interaction  of  insta¬ 
bility  (Tollmien-Schlichting)  waves.  The  instability  waves  are  due  to  the 
unstable  jet  stream.  It  is  therefore  desireable  to  study  the  interaction 
of  instability  waves  in  the  atmosphere  (due  to  the  existence  of  wind  shear) 
and  internal  gravity  waves  (due  to  atmospheric  stratification) .  In  order  to 
simplify  the  computation,  compressibility  will  be  neglected.  Since  the 
phenomena  under  study  do  not  depend  on  the  compressibility  in  any  crucial 
way,  this  seems  to  be  justified. 

Let  the  x-z  plane  be  horizontal,  the  y  axis  vertical  (positive 
upwards).  The  equations  of  motion  are: 

$  +7-(cv)  =  o,  o 

+(v-7)v]+  7 P  7Z$  (3 

if+(\/‘7)e-t c  •  (3 
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Here  p  is  the  density,  V  is  the  velocity,  p  is  the  pressure,  g  is  the  acceler¬ 
ation  of  gravity,  4  is  the  viscosity  coefficient,  0  is  the  temperature,  k  is 
coefficient  of  thermal  conductivity  and 

5  =  (0,  1,  0)4  (3.4) 

Koppel  (1964)  has  shown  that  both  viscosity  and  thermal  conductivity  must  be 
retained  in  order  that  the  small  disturbance  equations  be  non-singular  and 
thus  be  physically  meaningful.  It  can  be  shown  easily  that  these  equations 
are  equivalent  to 

&  •  =  Oj  (3.5) 

+  ( V’  7)c  ’  *  7  C  j  <3-6> 

f[$  +  (v-7)V]+7P+ffy^Vzv.  (3-7> 

In  deriving  these  equations  it  has  been  assumed  that  the  temperature  and 
density  are  related  by 

p/po  =  1  -  a  (0  -  eo),  (3.8) 

where  a  is  the  coefficient  of  thermal  expansion,  0q  is  a  reference  temperature 
and  pQ  is  the  density  at  temperature  0q;  the  Boussinesq  approximation  has  also 
been  made,  that  is,  in  the  momentum  conservation  equation  the  variation  of 
density  is  neglected  in  all  terms  except  the  one  containing  g. 
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Assume  that  the  flow  is  two-dimensional,  independent  of  z,  and 
introduce  a  stream  function  ¥,  with  u  and  v,  the  velocity  components  in  the 


x  and  y 

directions  given 

by 

um7F  • 

,  *■  — 

W 

Tx  • 

(3.9) 

Then  the  equati 

ons  of  motion  reduce  to 

(3.10) 

CfsF 

j  7 
<*X  7?  J 

(3.11) 

Now  the  linearized  equations  can  be  derived  by  writing 

£  $  +0(e  'l 

(3.12) 

r-r,  * 

* 

0(e')J 

(3.13) 

and  substituting  into  (3.10)  and  (3.11).  At  this  point  the  base  flow, 
given  by  pQ  and  ¥o,  is  taken  to  be  a  steady  flow  with  density  stratifica¬ 
tion  plus  a  wave  (either  a  stability  wave  or  an  internal  gravity  wave). 
Therefore  assume  that 


(3.14) 

(3.15) 
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The  undisturbed  density  is  given  by  POQ(y)  and  b'oo(y)  =  dYQO/dy  is  the  base 
shear  flow.  The  terms  Pq1  and  Y  ^  describe  a  wave,  of  amplitude  6,  propagat 
ing  through  the  flow.  In  order  to  have  a  consistent  description  one  must 
take 

</;  =  #.  (*.$*)  (*,&*)  *  O(S'). 

It  is  assumed  that 

Ij 

so  that  the  linearization  is  valid.  Substitute  into  the  basic  equations, 
again  using  the  Boussinesq  approximation.  The  quantities  poQ  and  Yoq 
satisfy  the  equations  identically.  For  the  others  it  is  found  that 


The  equations  satisfied  by  p1Q  and  Y10  are  identical  to  (3.19)  and  (3.20). 


(3.16) 


(3.17) 


(3.38) 


(3.19) 


(3.20) 
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If  wave  solutions  are  sought,  i.e.,  pQ1  and  or  p1Q  and  are 
taken  to  be 


a  - 


Uit) 


(3.24) 


(3.25) 


then  it  is  found  that 


i(-"  ■h*  Voo)  ^  -t-  ccip  (p  -  K  L.y  =  O 


(3.26) 


/  ('-w  *-*  Vtf,  )L  <p  -  <y  ^  K 


where 


^  tp-cWtf-yL'f.*; 


(3.27) 


_ 


££ 


(3.28) 


7 


and 


(3.29) 


Equations  (3.26)  and  (3.27)  govern  (in  the  linear  approximation)  the  propa¬ 
gation  of  waves  in  a  stratified  atmosphere  with  wind  shear.  In  various 
limits  some  well  known  results  are  obtained.  For  example,  ifg=g=K=0, 
there  is  a  homogeneous  atmosphere  in  the  absence  of  gravity  and  thermal  con¬ 
duction  and  (3.27)  becomes,  with  mr  Qfc 


(VO0-c)L  +(J$')  ll  Cp  =■  0J  (3.30) 
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which  is  just  the  Orr- Somme rfeld  equation  (Lin,  1955).  This  is  well  known 
in  the  theory  of  hydrodynamic  stability  and  describes  the  propagation  of 
stability  waves  in  shear  flows.  If,  on  the  other  hand,  Uqo  =  o  and  v  =  K  -  o, 
that  is,  an  inviscid  nonheat-conducting  fluid  is  considered,  with  no  wind, 
then  equations  (3.26)  and  (3.27)  reduce  to 

~yJ  <ps° .  (3-3i> 

This  is  just  the  equation  governing  the  propagation  of  internal  gravity  waves 
in  a  stratified,  incompressible  inviscid  fluid  (Tolstoy,  1963).  It  can  be 
noted  here  that  the  effect  of  the  Boussinesq  approximation  is  to  drop  a 
term  p  /4.  That  is,  the  exact  equation,  without  making  the  Boussinesq 
approximation,  would  be 

-JZZ  +  [ ~  if  f' ]<?=*.  (3.32) 

0 

If  the  stratification  is  small,  then  p  will  be  small  and  the  neglect 

2 

of  p  /4  is  reasonable. 

The  above  discussion  indicates  that  equations  (3.26)  and  (3.27) 
should  have  two  branches,  a  stability  wave  branch  and  an  internal  gravity 
wave  branch.  If  pQ1  and  YQ1  are  proportional  to  expCitOjX-w^t)]  and 
and  are  proportional  to  exp[i(a2X-W2t)]  then  equations  (3.21)  and  (3.22) 
show  that  p^  and  Y^  will  have  solutions  proportional  to  expCita^x-w^)) 
where 

°(j  =  Oi,*  ^  .  (3.33) 
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The  sketch  below  shows  schematically  the  dispersion  relation  for 
stability  waves  and  internal  gravity  waves.  It  is  clear  that  there  can  be  a 
resonance  between  various  waves.  For  example,  the  two  stability  waves  marked 
1  and  2,  and  the  internal  gravity  wave  marked  3,  will  satisfy  the  resonance 
conditions,  equation  (3.33). 

Stability 


In  order  to  proceed  any  further  with  the  calculation  a  model  must 
be  chosen  for  the  density  variation  and  velocity  field  in  the  undisturbed 
flow.  A  particularly  simple  model  has  been  chosen  so  as  to  minimize  the 
computation.  The  model  consists  of  a  duct  bounded  by  rigid  planes;  the 
Vcfiseflcf  frequency  is  assumed  to  be  constant;  that  is,  we  assume 

f  w 

is  constant;  and  finally  it  is  assumed  that  the  shear  flow  has  a  parabolic 
profile  (Uqo  =  1  -  y  )  since  this  is  the  simplest  profile  with  varying  shear. 
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In  all  that  follows,  dimensionless  variables  will  be  used.  The  height 
of  the  duct  is.^;  the  scaling  length  will  be^/2.  The  maximum  velocity  in  the 
duct  is  UQ  and  this  will  be  used  as  the  velocity  scale.  In  terms  of  these 
variables,  it  can  be  shown  that  the  linear  equations  (3.26)  and  (3.27)  are 

cr^  *  [-—  L  -tt*  -  e<K  <pj  <3-35> 

<rL  <?  =«[m-  vCtL]<p+±  L'cp  +uT^j  (3'36) 

where  a  is  the  complex  frequency, 

R  is  the  Reynolds  number,  =  UQf  /2)l 
P  is  the  Prandtl  number,  =  y/K 
J  is  the  Richardson  number,  =  egf/#20 

and  L  =  d2/dy2  -  a2. 

The  boundary  conditions  on  T]  and  ^  are  that 

%(*/)  =  (3-32) 

The  method  to  be  used  to  solve  equations  (3.35)  and  (3.36)  with 
the  boundary  conditions  (3.37)  is  to  choose  a  suitable  set  of  expansion 
functions  and  expand  1]  and  cp  in  a  series  of  these  functions.  There  is  no 
necessity  to  use  the  same  set  of  functions  to  expand  7]  as  is  used  for  Cp  and, 
in  fact,  different  sets  of  expansion  functions  will  be  used. 

A  suitable  set  of  expansion  functions  for  cp  is  the  set  of  eigen¬ 
functions  of 
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L-  LC„  i 

C„(±l)  =  (  ±/)-o. 


Since  equations  (3.38)  and  (3.39)  are  symmetric,  the  eigenfunctions  must  be 
either  symmetric  or  anti- symmetric .  For  the  present  calculation  it  is  only 
necessary  to  consider  the  symmetric  expansion  functions. 

The  symmetric  eigenfunctions  are 


where  y  is  the  n'th  root  of 
n 

tan  (i„)  -h<X  =  o  J 

+■  e(  lj 

and  c  is  a  normalization  constant, 
n 

This  approach,  expansion  in  a  series  of  orthogonal  functions, 
although  very  powerful  has  only  been  used  rarely  in  problems  of  this  sort. 
These  particular  expansion  functions,  the  Cn(y),  were  first  used  by  Dolph  and 
Lewis  (1958)  in  a  study  of  hydrodynamic  stability  and  were  later  suggested  by 
Eckhaus  (1965)  as  a  suitable  set  of  functions  to  use  in  non-linear  stability 
calculations. 


(3.38) 

(3.39) 


(3.40) 

(3.41) 

(3.42) 
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The  are  a  somewhat  unusual  set  of  functions.  Equations  (3.38) 

and  (3.39)  do  not  constitute  the  usual  Sturm-Liouville  (S.-L.)  problem.  The 
|Cnj  have  an  unusual  orthogonality  property,  which  is, 

“  Z.  ^  s  •  (3.43) 

Further,  the  do  not  have  the  usual  completeness  property  of 

the  S.-L.  functions,  (Courant-Hilbert,  1953).  The  S.-L.  functions  are  complete 
with  respect  to  the  space  of  functions  which  are  continuous  and  satisfy  the 
same  boundary  conditions.  If  a  function  does  not  satisfy  the  boundary  condi¬ 
tions,  then  the  expansion  coiveiges  non-uniformly  in  the  interior  and  differs 
from  the  function  by  a  constant  at  the  boundary. 

In  contrast,  it  is  easily  shown  that  if  one  attempts  to  expand 
cosh(oy)  or  sinh(ay)  in  a  series  of  the  |Cn|,  all  of  the  expansion  coefficients 
are  identically  zero.  This  implies  that  if  one  attempts  to  expand  an  arbitrary 
continuous  symmetric  function  in  terms  of  the  |c  \,  then  in  the  interior  the 
expansion  does  not  converge  to  the  function  but  to  the  function  minus  a  multiple 
of  cosh(ay) .  However,  it  can  be  shown  that  the  |cnJ  are  complete  with  respect 
to  the  space  of  continuous  functions  which  satisfy  the  boundary  conditions 
(3.37). 

A  suitable  set  of  expansion  functions  for  T|  are  the  eigenfunctions  of 


L  ^ 

G„  (-  t)  *  o . 


(3.44) 

(3.45) 


j. 

=  J  f(y)g(y)dy 

-l 
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Again  the  eigenfunctions  are  either  symmetric  or  anti- symmetric  and  only  the 
symmetric  functions  are  considered  for  the  present.  The  symmetric  eigen¬ 
functions  are 


(3.46) 


(3.47) 


Now  set 


/!*/ 


Af '/ 


(3.48) 


(3.49) 


and  take  the  inner  product  of  G^(y)  with  equation  (3.35)  and  C^(y)  with 
equation  (3.36).  After  integrating  by  parts,  using  (3.48)  and  (3.49)  and 
the  orthogonality  properties  equation  (3.43),  it  can  be  shown  that  equations 
(3.35)  and  (3.36)  reduce  to 

ra* =  ■  +**)  **  '•’**£.%,  v 

/t  n-  / 


50) 


ft  &  t 


(3.51) 


where,  with 


V,t  =  /- 


(3.52) 
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^  -  L y  y-2  cK  C„  \ 

'  J 

(fi,,  H S„,  Gm  X 


If  the  expansions  in  (3.48)  and  (3.49)  are  truncated  at  a  finite 
number  of  terms,  that  is,  if  T]  and  cp  are  approximated  by  a  finte  series,  then 
(3.50)  and  (3.51)  can  be  combined  to  give 


(3.53) 

(3.54) 

(3.55) 


(3.56) 

(3.57) 

(3.58) 

(3.59) 

(3. 60) 
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Then  a  is  just  an  eigenvalue  of  the  matrix  in  (3.56)  and  by  calculating  the 
eigenvalues  of  this  (2NX2N)  matrix  one  obtains  approximations  to  the  first  2N 
complex  frequencies.  Once  the  eigenvalues  are  known  the  wave  forms  can  be 
obtained  from  the  set  of  homogeneous  linear  equations  with  zero  secular  deter¬ 
minant. 

The  results  of  the  calculations  can  be  checked  in  a  number  of  special 
cases.  For  example,  if  J  =  0,  then  N  of  the  eigenvalues  ought  to  be  the  eigen¬ 
values  of  the  stability  modes  of  Poiseuille  flow.  This  has  been  verified  in  a 
number  of  cases. 

Numerical  results  are  presented  for  a  representative  case.  The 
values  assigned  to  the  various  scales  and  dimensionless  parameters  are: 

U  =  50  m/ sec . 
o 

}(  -  2  km. 

T  =^f/2U  =  20  sec. 

o  o 

4  2/  * 

v  =  5  x  10  cm  /sec. 

R  =  U  J?/2v  =  104 
o 

P  =  0.72 
J  =  0.005 

Since  (eq.  3.25)  the  stream  function  f  is  of  the  form 
»  =  <p(y>  e1(0K  +  <*>, 

then 

—  c $9+1  (o~ )  ■  OJj  (3.61) 

the  circular  frequency,  in  dimensionless  units  and 

*  The  kinematic  viscosity  is  taken  to  be  an  eddy  viscosity  for  which  this 
numerical  value  is  representative  of  atmospheric  conditions. 
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-  Jfti ?  (<r)  »  (3.62) 

the  damping  coefficient,  again  in  dimensionless  units. 

Thus, 

cp(y)  (3.63) 

In  Figure  3.1  the  variation  of  the  frequency,  u),  with  the  wave  number 
a  is  shown  and  in  Figure  3.2  the  variation  of  the  damping  coefficient,  6,  with 
wave  number  is  shown  (note  that  in  Figure  3.1  6  x  10  is  plotted).  Curves  are 

plotted  for  the  first  five  modes  where  the  modes  are  ordered  with  respect  to 
the  magnitude  of  the  damping  coefficient  for  small  0!.  It  is  obvious  that 
there  are  two  general  classes  cf  modes: 

1)  Modes  for  which  the  phase  velocity 

c  =  to/o  «  1,  (3.  64) 

these  may  be  called  the  ’’fast"  modes  and  include  modes  1,  3,  4  and  5. 

2)  Modes  for  which 

tt)  —^constant  as  Ct — >  »  (3.65) 

and  thus 

c  — ^  0  as  c t— ,  (3-66) 

these  may  be  called  the  "slow”  modes.  Mode  2  is  a  "slow"  mode,  and  although 
they  are  not  shown,  other  of  the  higher,  i.e.  more  rapidly  damped  modes  are 
also  "slow"  modes. 

As  discussed  above  internal-gravity  waves  can  propagate  in  a  station¬ 
ary,  stratified  fluid  and  Tollmien-Schlichting  or  instability  waves  can  propagate 
in  a  shear  flow  in  a  homogeneous  fluid.  Therefore  it  is  to  be  expected  that 
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there  will  be  both  internal-gravity  and  Tollinien-Schlichting  waves  in  a  strati¬ 
fied  shear  flow.  In  Table  3.1,  U)  and  6  have  been  listed  for  the  first  five 
modes  for  a  =  1.0,  and  also  cu  and  6  for  the  corresponding  case  in  the  homogen¬ 
eous  shear  flow  with  the  same  velocity  profile.  It  is  obvious  that  modes  1 
and  4  are  the  internal  gravity  modes  and  modes  2,  3  and  5  are  the  Tollmien- 
Schlichting  waves. 


TABLE  3 ■ 1 

FREQUENCY  AND  DAMPING  COEFFICIENT 
a  =  1.0,  R  =  104,  P  -  0.72,  J  =  0.005 


Stratified  Shear  Flow  Homogeneous  Shear  Flow 


!  Number 

(U 

6 

(D 

6 

1 

0.995 

0.0084 

2 

0.237 

-0.0039 

0.237 

-0.0032 

3 

0.964 

0.0358 

0.965 

0.0352 

4 

0.959 

0.0414 

- 

- 

5 

0.936 

0.0636 

0.936 

0.0633 
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It  can  also  be  shown  (see  eqs.  3.35  and  3.36)  as  R— and  assuming 

that  U  =  constant  =  1,  in  dimensionless  units,  that 
oo 

0)  =  a  -  0(J).  (3.67) 

This  shows  that  for  small  Richardson  number  (J)  the  internal  gravity  modes  are 
"fast11  modes  as  show'  in  Figure  3.1  and  Table  3.1. 

It  should  also  be  noted  that  for  0.8  <  a  <  1.1,  mode  2,  the  first  T-S 
mode,  is  unstable  as  indicated  by  the  fact  that  6  is  negative.  This  range  of 
wave  numbers  is,  to  within  the  accuracy  of  the  calculation,  the  same  range  of 
wave  numbers  for  which  the  homogenous  shear  flow  is  unstable.  The  Richardson 
number  is  so  small  that  the  stability  characteristics  are  not  appreciably 
affected.  For  larger  J  the  range  of  unstable  wave  numbers  decreases  and  for 
sufficiently  large  J  the  flow  is  stable  for  all  wave  numbers. 

Once  the  eigenvalues  were  obtained  it  was  a  simple  matter  to  calcu¬ 
late  the  wave  forms  from  the  set  of  homogenous  linear  equations  with  zero 
secular  determinant.  These  calculations  have  been  carried  for  a  number  of 
cases  and  some  sample  results  are  shown  in  Figures  3.3  -  3.6.  For  all  of 
these  cases  a  =  1.0.  In  Figure  3.3  the  nerturbation  velocity  field  is  shown 
for  mode  1  at  values  of  (oX  -  uut)  -  0,  (rr/4),  7tt/4  and  Figure  3.4  shows  the 

corresponding  results  for  mode  2.  The  velocity  scale  in  both  figures  is  the 
same.  The  same  information  is  presented  in  a  somewhat  different  way  in 
Figures  3.5  and  3.6.  Figure  3.5  shows  the  streamlines  for  mode  1  and  Figure 
3.6  the  streamlines  for  mode  2.  The  horizontal  scale  in  Figures  3.5  and  3.6 
has,  however,  been  compressed  in  order  to  show  a  wavelength  and  a  half,  that 
is,  0  <  (oX  -  uat)  <  3tt  in  Figures  3.5  and  3.6. 
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FIGURE  3.5  Streamlines  for  the  first  internal  gravity  wave  mode 

a  =  1.0,  R  =  104,  P  =  0.72,  J  =  0.005,  u>  =  0.995,  6  =  0.0084 


FIGURE 3.6  Streamlines  for  the  first  Tollmien-Schlichting  wave  mode 

a  =  1.0,  R  =  104  p  =  0.72,  J  =  0.005,  U>  =  0.237,  6  -  -0.0039 
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It  is  clear  that  the  internal  gravity  wave  is  largely  trapped  in  the 
high  velocity  region  at  the  center  of  the  jet  and  propagates  with  a  phase 
velocity  very  nearly  equal  to  the  maximum  velocity  of  this  shear  flow.  The 
T-S  wave,  on  the  over  hand,  is  spread  over  the  entire  duct  and  propagates 
with  a  phase  velocity  less  than  one-quarter  of  the  maximum  shear  flow  velocity. 

There  has  been  a  certain  amount  of  controversy  over  the  problem  of 
internal  gravity  wave  propagation  in  a  shear  flow.  The  crucial  point  is  a 
simple  one;  in  a  non-viscous,  non-heating  conducting  stratified  fluid  the 
first  order  perturbation  equations  are  singular  at  the  critical  point 
(the  point  where  the  phase  velocity  of  the  wave  equals  the  shear  flow  veloci¬ 
ty).  Although  Drazin  and  Howard  (1966)  in  their  comprehensive  review  of  the 
stability  of  parallel  flows  in  inviscid  fluid  were  careful  to  point  out  the 
existence  of  internal  gravity  waves  and  T-S  waves  in  stratified  shear  flows 
and  to  emphasize  that  the  relation  of  the  solutions  of  the  inviscid  problem 
to  the  solutions  of  the  viscous  problem  is  both  complicated  and  subtle,  other 
authors  have  not  been  as  careful.  One  difficulty  is  that  if  the  fluid  is 
stratified,  the  inclusion  of  viscosity  does  not  remove  the  singularity  as  it 
does  for  a  flow  in  a  homogenous  fluid.  However,  as  shown  by  Koppel  (1964), 
the  inclusion  of  a  non-zero  thermal  conductivity  does  remove  the  singularity. 
This  is  reasonable  since  the  density  stratification  is  due  to  a  temperature 
gradient.  The  equations  used  here  (3.35  and  3.36)  are  equivalent  (with  a 
different  notation)  to  Koppel's  final  result  (Koppel,  1964,  eq.  20).  It 
can  thus  be  seen  that  the  amplitude  of  the  motion  does  not  become  infinitely 
large  at  the  critical  layer.  In  the  neighborhood  of  the  critical  layer  the 
perturbation  velocity  does  reach  its  maximum  but  remains  finite. 
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The  major  results  are  that  this  stratified  shear  flow  supports  both 
internal  gravity  and  Tollmien-Schlichting  waves;  that  the  internal  gravity 
waves  are,  roughly,  trapped  within  the  layer  where  Uqo  >  c;  and  that  the 
phase  velocity  of  the  internal  gravity  modes  is  approximately  equal  to  the 
maximum  velocity  of  the  shear  flow. 
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4.  THEORETICAL  SIGNAL  PROCESSING 

4-1  Estimation  and  Detection  of  an  Unknown  Dispersive  Signal  Using  a 
Two-Dimensional  Array _ 

4.1.1  Introduction 

Consider  the  problem  of  detecting  and  estimating  an  unknown  tran¬ 
sient  signal  which  is  propagating  across  an  extremely  large  planar  array 
immersed  in  a  waveguide  such  as  the  ocean,  atmosphere,  or  the  earth's  surface. 
The  direction  of  the  propagation  and  the  phase  velocity  of  the  signal  are  un¬ 
known,  and  the  noise  field  is  assumed  to  be  Gaussian.  In  a  large  system 
tremendous  quantities  of  data  are  collected,  most  of  which  is  noise.  The 
transient  signals  have  a  relatively  long  time  between  arrivals  (such  as 
nuclear  events).  It  is  desirable  to  have  a  simple  procedure  for  detecting 
the  arrival  of  a  signal,  estimating  the  direction  and  phase  velocities,  and 
upgrading  the  data.  Once  the  data  has  been  segmented  into  periods  of  likely 
event,  the  optimal  methods  can  be  used  to  eliminate  false  alarms  and  improve 
the  signal  estimates. 

One  also  has  a  reception  problem  that  is  not  caused  by  noise.  The 
transmission  in  the  waveguide  is  dispersive  and  hence  the  phase  and  group 
velocities  depend  upon  frequency.  For  signals  having  a  wide  bandwidth,  it 
is  well  known  that  the  usual  methods  of  delay  and  sum  are  inadquate  for  array 
steering  since  the  signal  waveforms  are  systematically  different  at  different 
receivers,  Figure  4.1.  The  time  alignment  at  t^  is  incorrect  at  the  time  t ^ 
If  the  signal  is  modeled  as  a  section  of  a  stationary  stochastic 
process  in  time  and  space  dimensions,  the  problem  of  optimal  array  processing 
for  non-dispersive  signals  has  been  analyzed  by  several  investigators  (Faran 
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FIG.  4.1  Sketch  of  Signal  Arrivals  in  a  highly  dispersive  waveguide. 
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and  Hills,  1952;  Smith,  1956;  Berman  and  Clay,  1957;  Bryn,  1962;  Good,  1963; 

Burg,  1964;  and  Laster  and  Linville,  1966).  For  a  non-stochastic  non-dispersive 
signal  of  unknown  waveform,  direction,  and  phase  velocity,  a  type  of  least- 
squares  estimator  has  been  investigated  by  Levin  (1965).  If  the  phase  velocity, 
direction,  and  the  noise  covariance  function  are  known,  the  ordinary  least- 
squares  estimator  of  the  signal  is  maximum  likelihood  linear,  unbiased,  and 
minimum  variance  (Capon,  et  al.,  1967;  Graybill,  1961).  Levin's  estimator  is 
non-linear  since  it  is  found  by  searching  over  a  range  of  phase  velocities 
and  direction  cosines.  A  large  amount  of  computations  are  required  and  the 
optimal  array  steering  and  signal  processing  procedures  have  been  more  expen¬ 
sive  than  delay  and  sum  techniques. 

This  section  presents  a  relatively  simple  method  for  estimating 
the  phase  velocity  and  direction  of  an  unknown  plane  wave  signal.  The  finite 
Fourier  transform  is  applied  to  the  output  of  each  sensor  and  phases  of  the 
smoothed  frequency  components  are  calculated.  Ths  phases  are  linearly 
regressed  on  sensor  position  to  produce  estimates  of  the  wave  numbers.  By 
appropriate  smoothing  over  frequency  and  transforming,  estimates  are  obtained 
of  the  direction  and  phase  velocity.  The  precision  of  estimates  depends  on 
the  signal-to-noise  ratio,  the  square-root  of  the  number  of  sensors,  and  the 
number  of  wave  lengths  which  fall  in  the  array.  Since  the  procedure  involves 
Fourier  decomposition  of  the  signals,. it  is  easy  to  include  the  effect  of  dis¬ 
persion  across  the  array. 

4.1.2  Plane  Wave  Signal  and  Noise 

Assume  that  the  planar  array  has  K  sensors  which  are  located  at 
positions  (x^y^)  with  respect  to  the  coordinate  axis  x  and  y  (Figure  4.2). 
Let  Pk(t)  denote  the  sampled  output  from  the  k**1  sensor  for  t  =  0,  T,...,  (n-l)T. 
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Given  a  sampling  interval  t,  let  T  =  nT  denote  the  total  observation  period, 
which  is  the  same  for  all  sensors.  When  the  signal  is  present,  for  t=0,T, • • • ,T-t 


pR(t)  =  s  (t  - 


x^cos  0  +  yk  sin  0 


)  +  nk(t) 


(4.1) 


where  0  is  the  direction  of  propagation  and  v  is  the  phase  velocity  in 
the  x-y  plane.  The  noise  at  the  k^  sensor,  nk(t),  is  assumed  to  be  station¬ 
ary,  Gaussian,  zero  mean,  and  such  that  the  expected  value  E [nk(t)nk,(t) ]  =  0 
for  k  ^  k',  i.e.,  the  noise  is  incoherent. 

If  the  duration  of  the  signal  is  about  the  same  as  the  observation 
period  T  seconds,  then  the  signal  can  be  modeled  as  a  plane  wave  which  has 
the  following  finite  Fourier  series  representation  at  (x^>y^):  For  t-0 ,t, . . . ,T-t 


n-1 


s(t-Oxk-gyk)  =  £  S(u)^)eltuj^t_Q6ck  pyk^ 

d=o 

where  uk  =  ■—ft" ,  and  a(cu)  =  and  p(cu)  =  are  the  wave  slow¬ 

ness  components  in  the  x  and  y  directions.  For  a  given  j,  S(tu.)  is  the 

til 

complex  amplitude  of  the  uu .  frequency  component  of  the  signal  measured  at 

3 

(0,0),  the  origin  of  the  coordinate  system.  Since  s(t)  ir  real, 

S(u)  _.)  =  S*(<ju.),  where  the  star  denotes  complex  conjugate.  The  folding 
frequency  is  -|t  Hz. 

Let  Pk(<Uj)  be  the  j^1  Fourier  coefficient  of  Pk(t),  i.e., 


(4.2) 


for  n  =  T/t 

VV  =  7  l 


n-1 


mT 


m-0 
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From  (4.1)  and  (4.2)  It  follows  that  for  each  k  =  1,***,K  and  J  =  0,'"',n/2- 


pkCu»j)  _  |s(uij) |e3 Cujj(«xk+Pyk)  +  Y(ojj) 3  +  Nk^.j)  lpklelcf^ 

where  Y(u>.)  is  the  pliase  of  the  3*^  Fourier  component  of  the  signal  at  the 
J 

origin,  and  N^(ojj)  is  the  j  Fourier  coefficient  of  n^(t).  Letting 
2 

a  (id)  denote  the  power  spectrum  of  n  (t),  it  can  be  shown  that  N,  (u>)  has 
k  k  X 

a  mean  value  of  zero  and 


E(  |Nk(u>)  |2)  w|  02(<u) 

for  large  T,  i.e.,  the  j  coefficient  in  the  discrete  Fourier  repre¬ 
sentation  of  a  section  of  a  stationary  random  process  is  a  random  variable 
whose  mean  value  is  zero  and  whose  variance  is  the  area  under  the  power 
spectrum  in  a  narrow  band  of  bandwidth  i/T  Hz  centered  at  the  frequency 
j/T  Hz  (Hinich  and  Clay,  1968).  Moreover,  N,  (u>.)  and  N,  (m.,)  are 

K  j  K  j 

approximately  uncorrelated  for  j  4  j '  .  The  approximations  are  exact 

statements  if  the  process  is  white. 

The  complex  Fourier  coefficients  P,  (m.)  are  obtained  by  taking 

k  J 

the  finite  Fourier  transform  of  the  sections  of  the  sensor  outputs 
Pk(t),  u  <_  t  <  T.  The  fast  Fourier  transform  (Cooley,  Lewis  and  Welch, 
1967)  is  an  algorithm  which  provides  efficient  computer  calculation  of  the 
discrete  Fourier  coefficients.  These  complex  coefficients  will  contain  the 
pliase  characterirtics  of  the  signal. 


(4.3) 


(4.4) 


4-6 


Teledyne  Isotopes 


4.1.3  Signal  Processing 

Compute  the  P,  (to.)  for  frequencies  m.  such  that  S(tu.)  »  0. 
K  3  J  ,1 

For  each  sensor  k,  compute 


tan 


-1 


ImPk(uK) 

RePk(ujj) 


for  0  <  (jjj  <  tt/t,  provided  that  the  imaginary  part  of  P^ImP^,  is  smaller 
in  magnitude  than  RePk,  the  real  part  of  P^.  If  |ReP^|  <  | ImP^ | ,  then 
let 


TT 

2 


Repk(ouj) 


For  large  T  it  can  be  shown  by  Taylor  series  expansion  of  cp^ 
as  a  function  of  (Brownlee,  1965)  that  the  phase  of  P^Coj)  can  be 

written 


cpk(w)  =  a)[a(u))xk  +  p(cu)yk]  +  Y(u>)  +  ek(«>)  > 


where  ek(oo)  is  3  Gaussian  ’’•andom  variable  whose  expected  value  is 
E[ek(u>)]  =  0(T_1)  and  whose  variance  is 


<7k(oj) 

2T|S(U))|2 


3 


since  Repk  and  ImPk  are  uncorrelated. 


(4.5a) 


(4.5b) 


(4.6) 


(4.7) 
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Since  arctan  is  a  multiple  valued  function,  cp( uu)  is  defined  only 
up  to  an  added  factor  of  2nm.  One  can  determine  m  by  starting  at  low 
frequencies,  that  is  if  given  u>,  |(xk+1-xk)cos  6  (yk+]-yk)sin  0|  < 
where  \(oj)  is  the  wavelength  of  the  ujth  frequency  component,  then 

lvr\l  <n  +  lvrekl 

and  thus  there  is  no  problem  determining  the  proper  integer  multiple  of 
2n  for  each  k  provided  that  ak(<u)  TS  (uj)  is  small.  For  an  orthogonal 
array  it  is  sufficient  to  have  the  minimum  distance  between  sensors  to  be 
less  than  \/2.  For  short  wavelengths,  it  may  be  necessary  to  break  the 
array  into  smaller  subarrays  and  then  the  estimates  of  6  and  wave  slowness- 
derived  from  the  subarrays  can  be  averaged.  These  estimates  can  then  be  used 
to  estimate  the  proper  integer  multiple  of  2rr  for  each  sensor  phase.  This 
procedure  for  determining  the  m's  requires  more  processing  time  and  thus 
if  it  is  necessary  to  use  it,  the  least-squares  search  procedure  for  estimat¬ 
ing  direction  and  velocity  will  probably  be  at  least  as  efficient  as  the 
method  discussed  in  this  section.  However,  our  method  can  deal  with  disper¬ 
sion  whereas  the  least-squares  method  cannot  unless  the  dispe- sion  effect  is 
known  in  detail. 

let  cp(uj)  be  the  K  dimensional  vector  defined  by 
cp(tu)  '  =  [^(<1)),  •  •  •  ,cpK(<a)]  and  let  e(uy)  be  the  vector  e(u>)'  [e^uo) ,  *  * ' ,  eK(co)  ] 
where  the  prime  denotes  the  transpose.  Let  y(v')  be  the  vector  of  parameters 
y(u)) '-((uQf,ujg,Y)  and  let  A  be  the  K  x  3  matrix 
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[  xi  yi  A 


A  ^ 


•  • 


\.  XK 


J 


Writing  the  system  of  equations  defined  by  (4.6)  in  vector  and  matrix  form, 


cp(oj)  =  Ay(o))  +  e(<u)  . 


Letting  E  denote  the  K  x  K  diagonal  matrix  whose  kth  diagonal  element 
]  2 

is  "2'pa^(u,)j  the  weighted  least-squares  estimator  of  y  is 


A 


y(u>)  =  (A'E“JA)~]A’E  3cp(u>) 


where  ^(u>)'  =  (u$,u£,Y)  and  y  indicates  the  estimate  of  y  .  Since 


,-l 


E(c)  is  of  the  order  of  T  ,  it  can  be  shown  that 


A  /  1  \ 

E(y)  = 


The  linear  estimator  y  is  optimal  for  the  system  given  by  (4.9)  (Deutch, 
1965).  In  addition  to  the  estimates  of  wave  slowness  and  signal  phase, 
it  is  possible  to  obtain  an  unbiased  estimate  of  the  signal  energy 
| S ( cu)  |  •  From  (4.4)  and  the  fact  that  the  expected  values  of  ReN,.  and 
ImN^.  are  zero,  by  the  Jaw  of  large  numbers  it  follows  for  large  KT  that 

K 

i  L  lpkM2  *  | s(«>)  | 2  +  yt  ak(u,) 

k-1 


(4.8) 


(4.9) 


(4.10) 


(4.13) 


4-9 


Teledyne  Isotopes 


Thus 


K 


is(«,)i2  =!  £ci  pk(o,) 


7 


k=l 


2 

is  an  approximately  unbiased  and  consistent  estimator  of  | S( a>)  | 

A 

The  covariance  matrix  of  y  is  given  by 


(4.12) 


E[(y-y)(y-y)  ']  "  |S(u))|“2  (A f  S—1A)— 1 


(4.13) 


In  order  to  simplify  the  above  expressions,  let  us  restrict 
ourselves  to  the  case  when  the  noise  at  the  kth  sensor  is  white  over 
the  band  where  there  is  appreciable  signal  energy,  i.e.,  for  each 
k=l, • • • ,K 


a^(iu)  -  if  S(uj)  »  0 


(4.14) 


Further,  let  us  choose  the  origin  of  the  coordinate  system  such  that 
K  K 


k=l  ak 


Then  from  (4.8), 
(of  order  ) 


(4.10),  (4.11),  and  (4.14),  we  have  the  following  unbiased 
estimators  of  a,  (3,  and  Y  : 


(4.15) 
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a(io)  -  (ouD) 


-1 


r  K  K 

r  -2  2 

l  ak  yk 

L’k-1 


K  K 

1  \  \\|  -|  I  °k2Vk) 

k=l  /  '  k-1  / 


K 

1 

k-1 


(4.16a) 


P(u>)  =  (ujD) 


-1 


I  ffk2  *k)/  I  “k^AI 


’k=l 


k=l 


-  /  I  \  VkJ  | °v\%! 


(4.16b) 


where 


D  = 


(itf’di 

,  K 

(  T  - 2 

L  ak  yk 

'k-1  / 

k=l  j 

1  °k2  Vk 

k  l 


and 


A 

Y(cu) 


K  ,-1 


1 

l  k-1 


-2 


i<v 

k-1 


(4.17) 


Note  that  these  estimators  are  linear  functions  of  the  cp^  and  that  the 
linear  coefficients  do  not  depend  on  uu  . 

A  A 

The  following  expressions  for  the  variances  of  a  and  g 


and  the  covariances  can  be  derived  from  (4.13): 


A 

via) 


y 


2 

k 


k-1 _ 

2T<u2|SU)  |2d 


A 

v(e) 


K 


x 


2 

k 


k-1 _ 

2Tou2  |S(ou)  |2D 


(4.18a) 


(4.18b) 
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A  A 
COV  (a,e) 


l °k2  xk  4 

k-1 _ 

2Toj2  |S(c«)pD 


and 


A  A  A  A 

COV  (a,Y)  =  COV  (g,Y)  =  0  . 


(4.]8c) 


(4.18d) 


The  variance  of 


A 

V  (  y) 


A 

Y  is 


k-1 

2T|s(ou)  |2 


(4.19) 


In  order  to  better  understand  the  above  estimators  and  their 

2  2 

properties,  let  us  consider  an  example.  Let  \  ~  o  for  each  k  -  1, 

Let  the  array  have  a  cross  configuration  with  equal  spacing  of  distance  d 
between  the  sensors  in  each  wing  of  the  cross,  and  let  there  be  exactly  N 
sensors  in  each  arm  (Figure  4.3).  Assuming  there,  is  n£  sensor  in  the  middle, 
K  2N  and  the  distance  across  the  array  is  Nd .  Let  the  origin  of  the  co¬ 
ordinate  system  be  located  at  the  middle  and  let  the  arms  be  the  x  and  y 
axis.  Then 


K 


I  ’k'k  = 

k=l 


0 


and 


N/2 

2  ^j2d2=  N<l^p+2?..  d2  . 
J=1 


(4-20) 


(4.21) 
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FIGURE  4.3  Cross  array  with  sensors  located  on  coordinate  axis. 
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2  2 

Since  =  o  for  each  sensor,  (4.20)  implies  that  (4.15)  is  true,  and 
thus 


Thus  from  (4.16), 


o(u>) 


K 

k-1 

- R - 

\  2 

tt,Z*x  k 

k=l 


N/2 

12  l 

-N/2 


j  9  . 

Txj 


cudN  (N+l)  (N+2) 


A 

P(u>) 


K 

lw 

k-l 


"I* 

k-l 


N/2 

12  L  j 

.1=  -N/2 


uidN  (N+l)  (N+2) 


where 

and 

from 


cp  .  is  the  phase  measured  at  the  j**1  sensor  on  the 
is  the  phase  at  the  j*'1  sensor  on  the  y-axis  arm. 
(4.17),  the  estimator  of  the  signal  phase  function  is 


x-axis  arm 
Moreover 


k=l 


N/2 

1  V  ( 

2N  L  ^xj 

j=  -N/2 


(4.22a) 


(4.22b) 


(4.23) 
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From  (4.18)  and  (4.19)  the  variances  and  covariances  of  and 


A 

¥  are 


V(fi)  = 


T  u)  |S(o)>  |2  N  (N+l)  (N+2)  d' 


=  v(ft  , 


(4.24) 


, A  A.  .A  A.  .A  A 

cov  (fi,p)  =  cov  (a,Y)  =  cov  (g,y)  =  o 


(4.25) 


A 

V(Y)  = 


2KT  |S(«)|: 


(4.26) 


Approximating  N  (N+l)(N+2)  by  N  in  equation  (4.24),  and  since  the  wave- 
til 

length  of  the  u>  frequency  is 


\(u>)  =  -2=-^ 

O) 


x  =  x/cose 

x 


Xy  =  X/SIN0 


it  follows  that  the  proportional  standard  deviation  of  a  is 


sp(a)  _  (J  KT  tt  |s(ou)  I  x  cose  V 

a  \v  3  a  X  / 


(4.27) 


provided  that  0  /  ^  or  .  Thus  the  proportional  standard  deviation 

of  the  estimate  of  the  x-axis  component  of  wave  slowness  depends  inversely 
on  vt;  the  signal  to  noise  ratio  |S(iu)  |/a,  and  Nd/X  ,  the  number  of 
wavelengths  lying  within  the  length  of  the  x-axis  arm  of  the  array.  Clearly 
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the  proportional  standard  deviation  is  undefined  for  0  =  —  or  —  since 

A 

a  -  0.  A  similar  result  holds  for  SD(p)/p  with  \  sin  0  replacing 
\  cos  0.  We  have  shown  that  the  best  estimate  of  the  velocity  is  obtained 
from  an  array  arm  along  the  direction  of  propagation  of  the  signal. 


4.1.4  Estimation  of  Direction  and  Dispersion 

Let  n  denote  the  number  of  uu.  >  0 
s  3 

i.e.,  if  the  signal  is  narrow-band  then  n  /T  Hz 

s 

the  estimator  of  signal  direction  be  defined  by 


such  that  »  0 

is  the  bandwidth.  Let 


> 


A  _1 
0  =  tan 


(4.28) 


Where  the  sum  is  taken  over  the  oa.  which  defines  n  ,  if  the  term  in 

5  s’ 

the  bracket  is  less  than  one.  If  the  term  in  the  brackpt,  i.e.,  the  average 
A  a  A 

ratio  of  p  to  a  .  is  greater  than  one,  define  0  in  terms  of  the  arccot 


as  in  (4.5b). 

By  Taylor  series  expansion  of  (4.28),  it  follows  from  (4.11)  and 

/  N  A 

(4.18)  that  the  expected  value  of  0  is 


*  *  *  0  (kt) 


(4.28a) 


and  the  variance  is 


A 

V(0)  = 


(sin  0  cos  0) 


2 


2cov(a,e)  4 

oe  ,2 


(4.28b) 
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For  the  cross  array  example, 


A 

v(e)  = 


12a 


KTn2(Nd)2 


d^ 


4n2  |S(a>)  | 2 


If  the  signal  is  narrow-band  with  amplitude  S,  then  the  standard  deviation 
A 

of  0  is 


where  is  the  average  wavelength  of  the  band.  These  results  are  similar 

to  the  ones  obtained  by  Levin  for  the  least-squares  processor  (Levin,  1965). 

In  order  to  measure  the  signal. dispersion,  let  us  assume  that 
the  wave  slowness  6(u0  is  approximately  linear  in  u)  for  the  frequency 
range  of  interest,  i.e., 

6(u>)  =Vo2  +  p2  =  6o  +  fyo 


-2  2 


-2  2 


-2 


Farther  assume  that  £  a^  x^  =  jE  a^  yk  and  £  x^y^  =  0  ,  and  thus 
A  A  A  A 

from  (4.18),  V(a)  =  V(g)  ard  C0V(a,p)  =  0.  Defining 


A 

6(a)) 


it  then  follows  that  the  variance  of  6  is 


A 

V(6)  = 


-2  2 


2Td>  |  S(d>)  j  g  xk 


(4.30) 


(4.31) 


(4.32) 


(4.33) 
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SD(  V  _  US  _n§  Nd  r1 

60  ^  1  3  a  X0j 


(4.34) 


A  _2 

Thus  the  proportional  standard  deviation  of  6^  is  iuQ  times  the  proportional 

A 

standard  deviation  of  50  ,  which  implies  that  the  precision  of  the  estimate  of 
6^  is  less  for  low  frequency  narrow-band  signals. 


4.1.5  Signal  Detection 

A 

The  estimator  can  be  used  to  detect  the  presence  of  a  propa- 

A 

gating  signal.  Given  the  estimate  of  ,  the  95%  confidence  interval 

for  is 


A  A  A  A 

[fi  -  1.96  SD(fi0),  6q  +  1.96  SD(6q)3  (4-35) 

A 

where  SD(6q)  is  given  by  (4.34).  Since  6Q  >  0,  tlie  signal  is  detected 

A  A 

with  "95%  confidence"  if  6Q  >  1.96  SD(6q).  If  the  computation  required  for 

A 

the  estimation  of  6^  is  considered  too  lengthy,  the  estimator  6(u>q)  can 

be  used  to  detect  the  signal,  where  is  the  frequency  for  which  the  signal 

to  noise  ratio  is  a  maximum.  The  signal  is  detected  with  "95%  confidence"  if 

6  (u>0)  >  1.96  SD(fi)  where  SD(fi)  =  is  given  by  (4.33). 

A 

The  estimate  6^  of  6-^  ,  the  slope  of  6(®),  can  be  used  to 

test  whether  or  not  there  is  dispersion  in  the  signal.  The  95%  confidence 

A 

interval  for  6-j ,  given  6^,  is 

Cfij  -  1-96  SD(  6^ ),  +  1.96  SD^)]  , 

and  thus  if.  zero  lies  in  this  interval,  the  data  indicates  that  there  is  no 
dispersion. 
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4.2  Cophase:  An  Ad-Hoc  Array  Processor 


4.2.1  Introduction _ 

"Cophase”  is  an  ad  hoc  statistic  characterizing  signal  strength 
as  a  function  of  velocity  and  direction,  for  data  from  an  array  of  detectors 
receiving  propagating  signal(s)  and  uncorrelated  noise.  The  definition  of 
cophase  is: 


c(v,e)= 


N2  K  K— 1 

l  l  l 

n=N1  j=k+l  k=l 


Aj(“n>  +  Ak(“n)j  cos  [(»l  -  ‘k>  -  (‘j  -  *k>] 


w2  K  K-l 

l  l  l 

n=N  j=k+l 


Aj(“„>  4  Ak(“n>| 


where: 

V  is  the  assumed  phase  velocity; 

0  is  the  assumed  azimuth  from  which  signal  arrives; 

n  is  the  index  of  frequency  (frequencies  in  the  band  from  the 
Njth  to  the  N2th  harmonic  are  used); 

j  is  the  station  number  (there  are  K  stations) 

«Dn  is  the  n'th  frequency 

Aj(«i>n)  is  the  n’th  Fourier  amplitude  at  the  ,1'th  detector 

F 

ft.  is  the  n'th  Fourier  phase  at  the  j'th  detector 
J 

V 

ft,  is  the  phase  of  the  n'th  frequency,  arriving  at  the  j'th 
^  detector,  based  on  the  position  of  the  detector  and  the 
assumed  values  of  V  and  0 

Basically,  C(V,0)  calculates  the  normalized,  weighted  sum  over  all  —  K(K-l) 
detector  pairs,  and  over  all  frequencies,  of  the  cosines  of  the  difference 
between  the  Fourier  phase  differences  and  the  assumed  phases  differences. 
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Note  that  C(V,0)  is  a  statistic  with  an  expectation  of  0  if  the  K  records 
are  non-correlated .  It  has  an  expectation  of  +1.0;  if  the  K  records  represent 
a  non-dispersive  propagating  signal,  at  the  assumed  value  of  V  and  0  corre¬ 
sponding  to  the  actual  velocity  and  direction  of  the  signal.  Its  expectation 
at  other  assumed  (V,0)'s  would  approach  0  for  a  very-large-aperture  array  or 
for  a  broadband  signal,  but  may  in  general  have  non-zero  values  (side  lobes). 

In  applying  this  statistic  to  signal  detection  or  velocity/direction 
estimation,  the  first  step  will  be  a  Fourier  analysis  of  the  K  records;  then 
C(V,0)  can  be  calculated  for  each  point  in  a  grid  of  velocities  and  directions. 
If  the  value  of  C(V,0)  is  sufficiently  close  to  +1.0  for  some  (V,0),  then  a 
signal  has  been  detected.  The  point  (V,0)  at  which  C(V,0)  is  at  its  maximum 
is  an  estimate  of  the  corresponding  signal  parameters.  This  "cophase"  sta¬ 
tistical  method  of  detection  and  parameter  estimation  does  not  require  a 
knowledge  of  the  +  2nn  ambiguities  in  the  Fourier  phases,  or  of  the  origin 
time  of  the  signal  and  its  source  distance. 


4 . 2 .2  Distribution  under  the  Null  Hypothesis 

If  "sufficiently  close  to  +1.0"  is  to  be  used  as  a  criterion  for 
signal  detection,  the  word  "sufficiently"  will  require  a  quantitative  defini- 

p 

tion  based  on  the  statistical  distribution  of  C(V,0)  when  the  #  's  i»3'e 

« J 

uniformly-distributed,  independent,  random  variables  (i.e.,  when  the  K  records 

are  non-correlated  noise).  This  null  hypothesis  distribution,  calculated 

numerically,  is  discussed  below  for  the  case  in  which  A.(u)  )  =  1.0  for  all 

J  n 

j  and  n.  The  extension  of  these  results  to  A's  which  vary  with  n  should 

p 

be  trivial;  but  if  the  A's  as  veil  as  the  •  's  are  allowed  to  vary  randomly 
with  j,  the  results  will  require  some  re- interpretation. 
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Since  there  are  K  detectors  in  the  array,  there  will  be  K  inde- 

1 

pendent  values  of  $.(u>n)  for  each  frequency.  There  are  —  K(K-l)  different 

F  F 

values  of  the  term  (4.  -  4.)  ,  but  only  K-l  of  these  are  independent. 

3  K  ' 

p 

The  K-independent  4  •(%)  are  uniformly  distributed  from  0  to  2tt. 

1  F  F 

The  —  K(K.-l)  differences,  4j  -  4^  >  of  which  only  K-l  are  independent,  have 

a  triangular  probability  density  function;  its  maximum,  1/2tt,  is  at 

p  p  F  F 

#.-#.=  0  ,  and  it  decreases  linearly  to  0  at  4.-4.  =  +  2tt  •  From 

3  K  3  K 

F  F 

this  distribution  of  4 .  -  {.  ,  it  may  be  calculated  that  the  probability 

J  K 

F  F 

density  function  for  cos  (4  .  -  4. )  is 

3  K 


p  cos(4j  -  4^) 
•  * 


tt“^1-cos2(4j  -  4^) 


between  1,  and  0  elsewhere.  This  function’s  second  moment  about  its 

F  F  F  F 

mean  is  the  variance  of  cos  (4.  -  #.)  ,  where  4.  and  4  are  indepen- 

3  K  3  k 

dent  random  phases  uniformly  distributed  from  0  to  2tt;  the  variance  may 
be  shown  to  be  equal  to  0.5. 

If  C(q,V)  were  defined  in  such  a  way  that  the  two  inner  summations 
were  made  over  only  (K-l)  detector  pairs  with  (K-l)  independent  phase  differ¬ 
ences,  the  expectation  and  variance  of  C(9,V)  would  be  0  and  0.5/(K-l), 
respectively,  and  C(9,V)  would  be  distributed  approximately  normally  for  large 
K,  under  the  null  hypothesis. 

However,  C(9,V)  is  formed  by  a  summation  over  all  K  (K-l) 
detector  p'  vs.  If  the  ~  K(K-l)  phase  differences  were  independent  (but 
they  are  not),  C(9»V)  would  have  an  expectation  and  variance  of  0  and 
0.5,-|  K(K-l),  respectively,  and  C(9,V)  would  be  distributed  approximately 
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2  i 

normally  for  largo  K  ,  undpr  tho  null  hypothesis.  If  the  —  K(K-l)  phase 
differences  were  independent  (but  they  are  not),  C(@,V)  would  therefore  be 
a  much  more  oowerful  statistic  than  the  modified  C(0,V)  of  the  last  para¬ 
graph,  in  which  only  K-l  independent  phase  differences  are  used  at  each  fre¬ 
quency.  Since  only  K-i  of  the  ■—  K(K-l)  phase  differences  are  actually 
independent,  there  is  no  a  priori  reason  to  expect  C(g,V)  to  have  a  smaller 
variance  than  0.5/(K-l).  As  the  numerically  determined  distribution  of 
C(0,V)  will  show  however,  its  variance  is  not  only  smaller  than  0.5/(K-l); 
it  is  actually  equal  to  0.5/-iK(K-l) . 

Using  16,384  sets  of  10  random  phases  each,  the  distribution  of 
C(V,e)  was  computed  numerically  for  K  =  10  (10  stations)  and  N-j  =  N2 
(one  frequency  only),  assuming  that  each  detector  senses  only  uncorrelated 
noise.  This  null  hypothesis  distribution  is  shown  in  Figure  4.4  by  the 
curve  labelled  0.0.  This  distribution  has  a  variance  of  0.011  and  a  mean 
of  -0.001.  The  variance  is  exactly  equal  to  that  which  we  would  expect  if 
the  I  K(K-l)  terms  in  the  summation  were  independent  -  i.e.,  0.5/^(l0)  (10-1) . 
It  is  smaller  (and  the  test  is  therefore  more  powerful)  by  a  factor  of  5,  than 
if  we  had  used  only  K-l  independent  phase  difference  from  among  the  -|k(K-1)  . 
The  distribution  is  strongly  skewed;  its  cumulative  frequency  at  C(V,0)  =  0 
is  0.68,  its  maximum  is  between  -0.05  and  -0.10,  and  its  mean  is  -0.001 
(but  numerical  computation  errors  for  the  mean  may  be  as  large  as  0.001). 

4.2.3  Distribution  under  the  Alternative  Hypothesis 

The  distribution  was  re-evaluated  assuming  that  each  detector 
senses  signal  in  addition  to  the  uncorrelated  noise,  for  a  number  of  signal- 
to-noise  amplitude  ratios.  Several  of  these  are  shown  in  Figure  4.4  along 
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FIG.  4.4  Hie  probability  density  fraction  of  C(V,6),  evaluated  at  the  (V,e)  of 

the  sip&al,  using  a  single  frequency,  for  several  signal- to-noise  ratios. 
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with  the  null  hypothesis  curve  (signal-to-noise  ratio  of  0.0).  As  the  signal- to- 
noise  ratio  passes  in  excess  of  1.0,  there  is  a  rapid  transition  in  which  the 
distribution’s  variance  at  first  increases  and  then  approaches  zero,  as  its 
mean  approaches  1.0. 

The  curves  in  Figure  4.4  can  be  used  to  estimate  Ct  and  3,  the  proba¬ 
bilities  of  Type  I  errors  (false  alarm)  and  Type  II  errors  (undetected  signal), 
respectively,  and  to  determine  if  they  are  acceptably  small.  For  example, 
the  97.5  percentile  of  the  0.0  curve  and  the  75  percentile  of  the  1.0  curve 

are  both  0.5.  This  means  that  for  an  Q  of  2.5%  (acceptably  small),  3  would  be 

75%  (unacceptably  risky)  if  the  signal-to-noise  ratio  is  1.0;  but  3  decreases 
rapidly  to  0.2%  as  the  signal-to-noise  ratio  increases  only  slightly  to  1.414. 
Furthermore,  if  more  than  one  frequency  is  available  for  summation,  an  accept¬ 
able  value  of  3  is  possible  for  signal-to-noise  ratios  of  less  than  1.0.  For 
example,  if  16  frequencies  are  used  in  the  summation  ( N2~ N j  =  15),  an  a  of  2.5% 

would  correspond  to  a  3  of  1.0%  if  the  signal-to-noise  amplitude  were  only  .707. 

4.2.4  Evaluation  of  the  Cophase  Technique:  Numerical  Experiments 

To  evaluate  the  cophase  technique,  numerical  experiments  were  per¬ 
formed  using  synthesized  ’’signals”  and  "noise”  recorded  by  a  hypothetical  array. 
The  array  consisted  of  six  elements:  three  at  the  midpoints  of  the  358  km  legs 
of  an  equilateral  triangle,  and  three  at  points  one- fourth  of  the  distance  from 
each  vertex  to  the  next  (going  anti-clockwise).  This  configuration  was  studied 
by  Haubrich  (1968),  and  found  to  have  certain  optimum  characteristics. 

The  "signals”  were  synthesized  by  assuming  that  a  source  generated 
a  Dirac  delta- function  signal  at  a  distance  of  10,000  km  due  North  of  the 
array.  The  signal  propagated  dispersively  with  a  phase  velocity 

V(mps)  320  -  0.5  sinh  (  1Mt  -  232  ) 
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The  corresponding  group  velocity  exhibits  a  maximum  at  f  =  232/736  cycles/min, 
or  a  period  of  slightly  more  than  three  minutes.  The  resultant  synthetic 
signals  (Fig.  4.5)  are  reasonable  facsimiles  of  actual  microbarograph  records 
of  acoustic-gravity  waves  from  nuclear  explosions  (Dorm  and  Shaw,  1967). 

The  4-hour  synthetic  signals  at  the  six  hypothetical  array  points 
were  next  Fourier  analyzed;  then  their  cophase  function  was  calculated  using 
various  parts  of  the  signals'  spectra,  for  velocities  between  250  and  400 
mps,  and  for  360°  of  azimuth. 

dsing  only  N  =  97  (ninety-seventh  harmonic  -  2.4  minute  period) 
in  the  frequency  summation,  the  cophase  (Fig.  4.6)  has  a  peak  at  (V,0)  (320,0°), 

corresponding  to  the  actual  signal.  The  peak's  amplitude  is  0.994  -  almost  per¬ 
fect  cophase  (but  not  quite;  since  the  signal  is  dispersive).  However,  there 
are  several  closely- spaced  side  lobes,  the  highest  of  which  are  0.976  at 
(V,0)  =  (310, ±100°).  Since  the  wavelength  at  this  frequency  is  one-half  the 
detector  spacing,  the  2nn  phase  errors  in  the  Fourier  analysis  produces  closely 
spaced  side  lobes  in  the  cophase  function.  This  phenomenon  may  also  be  regarded 
as  a  ''ringing”  in  velocity-direction,  due  to  the  extreme  narrowness  of  the 
filter  (only  one  frequency  was  used).  By  summing  over  frequency  from  N  -  95 
to  N  =  100,  the  largest  side  lobe  drops  in  amplitude  only  to  0.868,  since  the 
frequency  filter  is  still  extremely  narrow.  Thus,  summing  over  adjacent  fre¬ 
quencies  within  a  very  narrow  band  is  not  useful  in  the  noise-l'ree  case. 

Using  only  N  =  33  (7.2  minute  period)  in  the  summation  (Fig.  4.7) 
the  side  lobe  reaches  a  height  of  0.997  -  the  same  as  the  main  peak  -  at 
(V,0)  =  (400, ±180°).  The  largo  spacing  of  the  lobe  from  the  main  peak  is 
consistent  with  the  much  larger  wavelength. 
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Cophase  of  the  7.2  min.  component  of  the  synthetic  signal  detected 
by  the  6-element  array. 
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Using  N  =  30  through  N  -  100  (Fig.  4.8),  the  main  peak  of  0.993 
clearly  dominates  over  the  region.  The  highest  side  lobes  are  only  0.133, 
at  (V,0)  =  350,  _+  50°).  The  broad  band  of  frequencies  integrated  provides 
effective  cancellation  of  the  side  lobes  appearing  at  different  positions  in 
the  cophase  functions  of  individual  frequencies.  The  result  of  using  a 
wide  band  of  frequencies  is  the  high  resolution  of  the  high  frequencies 
(Fig.  4.8)  combined  with  even  less  ambiguity  than  the  low  frequencies  alone 
(Fig-  4.7). 

When  computer-generated  noise  is  added  to  the  set  of  six  synthetic 
signals,  the  cophase  function  begins  to  degenerate.  Figure  4.9  shows  the 
cophase  calculated  with  one  frequency  (N  =  50),  and  a  signal- to-noise  ampli¬ 
tude  ratio  of  about  0.7.  There  is  no  peak  at  (V,0)  =  320,  0°),  where  the 
signal  should  be.  There  are  seven  peaks  with  cophase  exceeding  0.4,  the 
highest  of  which  is  0.548.  These  represent  random  noise  reinforcement,  and 
coincide  with  neither  the  signa~  nor  its  side  lobes. 

However,  if  cophase  is  recalculated  summing  over  frequency  from 
N  =  30  through  N  =  100,  a  peak  of  0.402  appears  at  (V,{))  320,  0°),  corre¬ 

sponding  exactly  to  the  signal  (Fig.  4.10).  In  this  case,  however,  the  noise 
is  sufficiently  cancelled  by  the  summation  over  71  frequencies,  so  the  highest 
peaks  beside  the  signal  itself  has  a  cophase  of  only  0.106  at  ( V, 0)  =  360,  -40°). 
The  result  may  be  unambiguously  interpreted  as  representing  a  signal  from  due 
north  with  a  phase  velocity  of  320  mps. 

A  final  pair  of  numerical  experiments  were  performed  to  evaluate 
the  use  of  the  cophase  method  in  conjunction  with  the  actual  Teledyne-Isotopes 
array.  First,  the  same  signal  described  above  was  synthesized  for  each  point 
in  the  array.  A  two-hour  sample  was  analyzed  using  N  =  15  (8  minute  period) 
t irough  N  =  50  (2.4  minute  period).  The  result  (Fig.  4.11)  shows  a  well- 
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FIG.  4.8  Cophaae  of  the  2.4  -  8.0  min.  band  of  the  synthetic  signal  detected 
by  tiie  6-element  array. 
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Cophase  of  the  4.8  min.  component  of  mixed  synthetic  signal  and  noise 
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defined  peak  of  cophase  0.998  at  N,8  -  (320,0°).  There  are  no  side  lobes  ex¬ 
ceeding  0.10  cophase  in  the  range  of  velocities  from  200  -  800  mps.  , 

The  eecond  of  this  last  pair  of  numerical  experiments  began  with 
the  synthesis  of  signals  similar  to  those  used  above,  but  with  a  mean  phase 
velocity  of  640  mps  at  a  period  of  12  minutes.  The  analysis  of  the  resulting 
two-hour  record,  based  on  frequencies  from  N  -  8  ( 15-minute  period)  to  N  -  12 
(10-minute  period),  is  shown  in  Figure  4.12.  Because  the  wavelengths  are 
larger  than  the  entire  array,  and  only  six  frequencies  are  used  in  the  summa¬ 
tion,  the  1.00  peak,  which  corresponds  to  the  signal,  affords  poor  velocity 
resolution,  but  the  half-power  points  are  at  less  than  ±10°  of  azimuth  from 
the  main  peak.  The  relatively  poor  velocily  resolution  is  an  intrinsic  proper¬ 
ty  of  the  array  at  these  large  wavelengths. 

4.2.5  Conclusions 

The  cophase  method  is  a  fast  method  of  detecting  signals  of  down 
to  one-half  the  noise  amplitude,  and  of  estimating  the  velocity  and  direction. 
The  most  important  limitations  are  the  intrinsic  limitations  of  the  array  size 
and  pattern,  and  of  the  signal  duration  and  bandwidth  but  not.  of  this  statis¬ 
tical  procedure,  pe r  se . 

The  method  is  presently  used  routinely  in  the  analysis  of  Teledyne 
Isotopes  microbarograph  and  doppler  array  data-  Several  interesting  prelimin¬ 
ary  results  not  obtainable  by  older,  simpler,  more  intuitive  techniques,  have 
already  been  arrived  at  by  the  use  of  the  cophase  method.  These  will  appear 
in  future  project  reports. 
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5.  SUMMARY  AND  CONCLUSIONS 

The  accomplishments  of  the  research  effort  during  the  period  1  September 
1968  to  30  September  1969  can  be  summarized  as  follows. 

1.  A  program  to  adjust  and  recalibrate  all  the  microbarographs  and  obtain 
their  phase  response  and  coherence  with  respect  to  a  ’’standard"  instru¬ 
ment  was  carried  out  when  the  system  became  operational. 

2.  Three  new  phase-path  sounder  systems  were  built,  which  added  to  the 
one  in  existence  made  it  possible  to  set  up  a  4-element  array  which, 
although  rudimentary,  would  permit  phase  velocity  and  direction  esti¬ 
mates  of  ionospheric  disturbances  to  be  obtained.  The  sounding  frequency 
used  in  this  doppler  array  was  4.824  MHz;  the  average  reflection  heights 
corresponding  to  this  frequency  are  about  170  km  during  daytime  and  200 
km  during  nighttime.  The  doppler  instruments  were  built  after  Davies 
(1962)  with  some  modifications  introduced  to  meet  our  specifications 

and  mode  of  operation.  Since  automatic  digitizing  equipment  was 
already  available  at  four  stations  of  the  microbarograph  array, 

(Catskill,  Thomhurst,  Lebanon  and  Westwood),  it  was  decided  to  locate 
the  doppler  receivers  at  these  stations  and  the  transmitter  at  Sterling 
Forest,  N.Y.  This  configuration  gave  an  average  separation  between 
ionospheric  reflection  points  of  about  80  kilometers.  The  data  from 
the  doppler  array  were  fed  into  spare  channels  of  the  data  acquisition 
system  and  recorded  together  with  the  geomagnetic  and  pressure  data  on 
incremental  digital  magnetic  tape. 
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3.  Numerous  calculations  were  undertaken  to  compare  the  theoretical 
response  of  the  Teledyne  Isotopes’  microbarograph  array  with  those 
of  a  number  of  other  arrays;  a  10-element  isometric  array,  a 
15-element  isometric  array,  a  21-element  isometric  array  and  the 
LAMA  array.  These  calculations  showed  that  for  wavelengths 
between  60  end  300  kilometers  the  response  of  the  Teledyne  Isotopes’ 
array  (10  elements)  is  very  similar  to  that  of  the  IAMA  array 

(20  elements)  and  that  of  a  21-element  isometric  pi  ray.  It  there¬ 
fore  follows  that  not  much  improvement  in  array  response  could  be 
expected  by  adding  new  instruments  or  changing  the  spatial  arrange¬ 
ment  of  the  present  array. 

4.  Atmospheric  disturbances  generated  by  the  Chinese  thermonuclear 
explosion  of  29  September  1969  were  detected  by  the  microbarograph 
array.  These  effects  showed  up  in  two  frequency  ranges  corresponding 
to  periods  of  2-5  minutes  and  12-18  minutes  respectively.  In  the 
shorter  period  range  the  arrival  corresponding  to  the  short  groat- 
circle  path  was  recorded,  while  in  the  longer  period  range  the 
arrival  corresponding  to  the  long  great-circle  path  was  recorded. 

The  group  velocities  of  these  arrivals  were  of  the  order  of  300  m  sec 
and  550  m  sec  *  respectively.  Some  conclusions  can  be  drawn  from  the 
comparison  of  these  observations  and  those  of  Tolstoy  (1967,  1968). 
a)  short  period  waves.  These  waves  arrive  with  group  velocities 
of  the  order  of  300  m  sec”*.  The  center  frequency  of  the  events 
compared  takes  different  values  in  each  case  ranging  from  0.002  Hz 
(8  minutes)  to  0.0067  Hz  (2.5  minutes);  that  is,  the  bandwidth  of 
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these  arrivals  is  of  the  order  of  0.00167  Hz  -  0.0167  Hz 
(2  minutes  -  10  minutes).  The  average  phase  velocities  are 
300  ±  10  m  sec-*  and  the  wavelengths  range  from  36  km  to 
180  km.  The  arrival  corresponding  to  the  long  great-circle 
path  cannot  in  general  he  recognized. 

b)  long  period  waves.  The  group  velocities  of  these  waves  are 
of  the  order  of  550  m  sec  *  and  the  phase  velocities  of  the 
order  of  570  ±  30  m  sec-*.  The  center  frequencies  also  vary 
from  event  to  event  as  in  the  short  period  case.  For  the 
events  compared  the  center  frequencies  were:  3.64  CPH  (16.5 
minutes),  4.62  CPH  (13  minutes),  and  4.0  CPH  ( 15  minutes) 
respectively,  that  is,  the  bandwidth  extends  from  0.000925  Hz 
to  0.00138  Hz  (18  minutes  -  12  minutes).  In  general,  both  the 
short  and  the  long  great-circle  path  arrivals  can  be  identified, 
however,  the  short  great-circle  path  arrival  corresponding  to 
the  29  September  1969  explosion  could  not  be  extracted  from  the 
background  while  the  long  path  arrival  was  identified  with  very 
little  filtering. 

c)  It  is  probable  that  results  in  a)  and  b)  above  are  an  indication 
of  the  effects  of  the  seasonal  variation  of  the  neutral- temperature 
structure  on  the  propagation  of  these  waves  since  the  events  com¬ 
pared  took  place  in  different  seasons,  i.e.,  during  late  spring, 
mid-summer  and  fall. 
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5.  Analyses  of  the  phtse-path  doppler  records  corresponding  to  the  French 
test  of  24  August  1968  showed  that  at  the  time  at  which  the  600  m  sec  * 
pressure  sigoal  was  identified,  there  was  a  group  of  doppler  oscilla¬ 
tions  with  an  amplitude  of  1  Hz  and  periods  of  10  -  15  minutes.  The 
signal- to-noise  ratio  was  3:1.  Since  only  one  doppler  sounder  was  in 
operation  at  that  time,  no  information  could  be  obtained  concerning 
the  speed  and  direction  of  travel  of  this  signal.  This  finding  was 

a  very  important  factor  in  redirecting  the  emphasis  of  the  project 
towards  ionospheric  studies.  Unfortunately,  during  the  time  corre¬ 
sponding  to  the  arrival  of  infrasonic  signals  in  the  September  1969 
test  a  magnetic  storm  was  in  progress  and  the  doppler  background  level 
was  very  high.  It  is  possible  that  the  doppler  signals  have  been 
masked  by  the  large  phase-path  fluctuations  induced  by  the  storm. 

6.  Observations  of  phase-path  variations  during  the  Satum-Apollo  IX 
launch  suggest  that  some  acoustic  energy  in  the  period  range  0.5  -  5.0 
minutes,  generated  by  the  rocket  engine,  is  channeled  at  ionospheric 
levels  as  evidenced  by  the  phase-path  doppler  records.  Analyses  of 
microbarograph  records  for  this  event  did  not  show  any  signal  in  this 
period  range. 

7.  A  preliminary  study  of  ionospheric  background  motions  was  carried  out 
using  doppler  and  magnetometer  records.  Power  and  cross  spectral 
analyses  of  these  records  indicated  that  the  background  activity  could 
be  broadly  separated  into  two  period  ranges:  periods  longer  and  shorter 
than  5  minutes.  The  shorter  periods  correlate  often  with  magnetic 
micropulsations.  The  analyses  of  these  events  showed  high  coherence 
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between  the  NS  magnetic  component  and  the  phase-path  variations,  sug¬ 
gesting  that  the  ionospheric  motions  producing  the  phase-path  changes 
are  due  to  the  modulation  of  the  ionospheric  plasma  by  transverse 
hydromagnetic  waves  *  Analyses  of  microbarograph  records  correspond¬ 
ing  to  these  time  intervals  of  high  correlation  did  not  show  any 
unusual  infrasonic  activity.  In  the  longer  period  range  the  power 
spectra  of  phase-path  sounder  records  obtained  simultaneously  at  three 
sounding  frequencies  (corresponding  to  three  ionospheric  heights) 
showed  systematic  peaks  and  slope  breaks  which  appear  to  shift  towards 
longer  periods  for  greater  reflection  heights  (higher  sounding  fre¬ 
quencies).  The  frequencies  corresponding  to  these  spectral  peaks  had 
values  of  the  same  order  of  the  local  value  of  the  Vb'istilA'  frequency 
for  an  adiabatic  atmosphere.  These  findings  lead  to  the  conclusion  that 
a  large  portion  of  the  ionospheric  motion  background  for  periods  shorter 
than  5  minutes  has  a  hydromagnetic  origin  while  periods  longer  than  5 
minutes  are  probably  due  to  internal  gravity  wave  activity. 

8.  Investigation  of  the  spatial  coherency  of  the  signal  and  noise  fields 
for  the  array  of  microbarographs  was  carried  out  using  the  signals  of 
the  French  test  of  24  August  1968  (2-5  minute  period  range)  and  compar¬ 
able  data  samples  of  background.  The  results  indicated  that: 

a)  Acoustic  signals  have  coherencies  well  above  the  background  noise 
at  separations  greater  than  3  or  4  kilometers. 

b)  The  signal  coherence  remains  high  out  to  separations  of  50  kilometers. 

c)  Transverse  signal  coherencies  are  higher  than  parallel  signal  coher¬ 
encies,  that  is,  the  waves  are  long-crested. 
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9.  A  theoretical  study  of  the  interaction  of  instability  waves  in  the  at¬ 
mosphere  (due  to  existence  of  wind  shear)  and  internal  gravity  waves 
(due  to  atmospheric  strati f teat  ton)  was  carried  out  by  Dr.  C.E.  Grosch. 
The  model,  whtch  neglects  compressibility,  consists  of  a  du't  bounded 
by  rigid  planes,  the  Vtfis&lg  frequency  is  assumed  to  be  constant,  and 
the  shear  flow  has  a  parabolic  profile,  |uoo  =  1  -  y2j  .  The  major 
results  are  that  this  stratified  shear  flow  supports  both  internal 
gravity  and  instability  waves;  that  the  internal  gravity  waves  are, 
roughly,  trapped  within  the  layer  where  UQo  >  c  ;  and  that  the 
phase  velocity  of  the  internal  gravity  modes  is  approximately  equal 
to  the  maximum  velocity  of  the  shear  flow. 

10.  In  the  theoretical  aspect  of  signal  processing  two  methods  have  been 
developed . 

a)  The  first  one  deals  with  the  estimation  and  detection  of  an 
unknown  dispersive  signal  using  a  two-dimensional  array  and 
represents  the  effort  by  a  member  of  our  group  (Dr.  M.  Hinich), 
to  find  a  simple  method  for  estimating  the  phase  velocity  and 
direction  of  a  plane  wave  signal,  buried  in  noise,  propagating 
across  a  two-dimensional  array.  The  method  uses  a  linear  re¬ 
gression  on  sensor  position,  of  the  Fourier  phases  of  the  sensor 
outputs,  in  order  to  produce  estimates  of  the  wave  number.  By 
appropriate  smoothing  over  frequency,  estimates  of  the  direction 
and  phase  velocity  are  obtained.  The  precision  of  the  estimates 
depends  on  the  signal-to-noise  ratio,  the  square  root  of  the 
number  of  sensors,  end  the  number  of  wavelengths  which  fall  in 
the  array. 
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b)  The  second  method  represents  the  effort  of  another  member  of  the 
group  (Dr.  E.  Posmentier)  and  consists  of  an  ad  hoc  array  processor 
that  characterizes  signal  strength  as  a  function  of  phase  velocity 
and  direction.  This  statistical  procedure  which  is  called  Cophase 
Analysis  throughout  this  work,  basically  calculates  the  normalized, 
weighted  sum  over  all  1/2  K  (K-l)  detector  pairs,  and  over  all 
frequencies,  of  the  cosines  of  the  difference  between  the  Fourier 
phase  differences  and  the  assumed  phase  differences.  Several  tests 
of  this  program  with  synthetic  and  real  data  have  shown  that  the 
cophase  method  is  a  fast  method  of  detecting  signals  of  down  to 
one-half  the  noise  amplitude,  and  of  estimating  the  velocity  and 
direction.  The  most  important  limitations  are  the  intrinsic 
limitations  of  the  array  size  and  pattern,  and  of  the  signal  dura¬ 
tion  and  bandwidth,  but  not  of  this  statistical  procedure  per  se. 

To  conclude,  the  authors  wish  to  recommend  that,  since  the  acoustic  energy 
generated  by  nuclear  explosions  in  the  atmosphere  is  distributed  with  height  accord¬ 
ing  to  various  propagating  modes,  observations  with  both  microbarograph  and  pi ase- 
path  sounders  be  continued  and  that  an  effort  be  made  to  expand  the  existing  phase- 
path  sounder  array  to  operate  at  two  sounding  frequencies. 
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A  suimary  is  given  of  significant  results  obtained  during  the  research  effort 
corresponding  to  the  period  1  September  1968  -  30  September  1969.  Particular 
emphasis  was  given  to  the  construction  of  a  phase-path  doppler  sounder  array. 
Signals  detected  by  the  microbarograph  array  on  the  occasion  of  the  Chinese  nuclear 
test  of  29  September  1969  are  compared  with  signals  from  previous  tests.  There 
appears  to  be  a  seasonal  ef  fect  on  the  propagation  and/or  bandwidth  of  these  sig¬ 
nals.  Preliminary  work  on  ionospheric  motion  background  using  phase-path  sounder 
data  indicates  that  the  background  activity  can  be  broadly  separated  into  two 
period  ranges:  periods  longer  and  shorter  than  5  minutes.  The  shorter  periods 
appear  to  have  a  hydromagnctic  origin  while  the  longer  periods  are  probably  due 
to  internal  gravity  wave  activity.  A  study  of  the  spatial  coherency  of  the  signal 
and  noise  field  for  the  array  of  mici obuiugrapha  ir.dicotcd  that  sc  uuu .. c  c i gnul » 
have  coherencies  well  above  the  background  noise  at  separations  greater  than  ‘i  or 
4  km.  The  signal  coherency  remains  high  up  to  separations  of  50  km  and  transverse 
signal  coherencies  are  larger  than  parallel  signal  coherencies  as  would  be  expected 
from  long-crested  waves.  A  theoretical  investigation  of  the  interaction  between 
instability  waves  and  internal  giivity  waves  in  the  atmosphere  showed  that  internal 
gravity  waves  are  trapped  within  a  layer  where  the  shear  flow  velocity  is  greater 
than  the  speed  of  sound  and  that  the  phase  velocity  of  the  internal  gravily  waves 
is  approximately  equal  to  the  maximum  velocity  of  the  shear  flew.  As  part  of  .the 
theoretical  aspect  of  signal  processing  two  array  processing  techniques  have  been 
developed . 
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